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ABSTRACT 


Although  ring-on-ring  test  results  have  been  used  in  the  past  to  obtain 
information  on  the  strength  of  glass,  no  methodology  has  so  far  been  developed 
in  the  literature  explicitly  relating  such  results  to  the  load  capacity  of 
cladding  glass.     The  main  purpose  of  this  report  is  to  propose  such  a  method- 
ology.    The  proposed  methodology  makes  use  of  recent  advances  in  the  modeling 
of  the  fracture  mechanics  behavior  of  glass  and  the  calculation  of  stresses  in 
plates  exhibiting  geometric  nonlinearity .     Evidence  is  presented  which  strongly 
suggests  that  the  probability  distribution  of  the  load  capacity  of  cladding 
glass  panels  whose  failure  is  due  to  surface  flaws  can  be  estimated  reliably 
on  the  basis  of  results  of  ring-on-ring  tests  used  in  conjunction  with  (a) 
numerical  methods  for  the  analysis  of  stresses  in  plates,  and  (b)  information 
on  the  elastic  and  fracture  mechanics  behavior  of  glass  currently  available  or 
that  can  be  obtained  routinely.     Two  interesting  findings  are  noted.  First, 
owing  to  the  way  in  which  results  of  ring-on-ring  tests  are  utilized,  the 
relatively  large  variabilities  typical  of  fracture  mechanics  parameters,  as 
well  as  the  uncertainties  with  respect  to  the  shapes  of  surface  flaws,  have  a 
minor  effect  on  the  estimation  of  load  capacities.     Second,  two-parameter 
Weibull  distributions,  previously  used  in  the  literature  to  model  the  strength 
of  glass  and  the  load  capacity  of  cladding  panels,  are  not  consistent  with 
experimental  results.     On  the  other  hand,  three-parameter  Weibull  distributions 
model  the  observed  glass  behavior  credibly. 

Keywords:     buildings;  engineering  mechanics;  failure;  fracture  mechanics; 

glass;  loads  (forces);  probability  theory;  ring-on-ring  tests; 
strength. 
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1.  INTRODUCTION 


For  most  common  construction  materials  (e.g.,  steel  or  concrete),  the  load 
capacity  of  structural  members  is  determined  on  the  basis  of  strength  data 
obtained  from  tests  conducted  on  small  size  standard  specimens,  rather  than 
from  the  destructive  testing  of  full-size  members.     The  economy  inherent  in  the 
use  of  small-size  standard  specimens  is  due  not  only  to  their  lower  cost  as 
compared  to  the  cost  of  full-size  members,  but  also  to  the  fact  that  strength 
data  obtained  by  testing  a  sufficiently  large  number  of  such  specimens  can  be 
used  for  designing  a  wide  variety  of  structural  members  with  different 
configurations,  types  of  loading,  and  sizes. 

Procedures  for  determining  the  load  capacity  of  glass  panels  on  the  basis  of 
strength  data  obtained  by  testing  small  standard  specimens  do  not  currently 
exist.     For  this  reason  design  charts  issued  by  glass  manufacturers  have  tradi- 
tionally been  based  on  destructive  tests  performed  on  full-size  glass  panels 
(see,  e.g.,  reference  1).     It  has  been  pointed  out  in  the  literature  (see,  e.g., 
reference  2),  that  such  charts  exhibit  significant  inconsistencies.  Such 
inconsistencies  are  due  at  least  in  part  to  the  relatively  small  numbers  of 
panels  (between  2  and  30 — see  references  4,  5,  and  6)  used  in  most  of  the 
tests  on  which  the  charts  were  based.     Since  the  charts  cover  a  wide  range  of 
panel  sizes,  the  number  of  panels  that  would  have  to  be  subjected  to  destructive 
tests  in  order  to  develop  dependable  data  for  design  could  in  practice  be 
prohibitive. 

It  therefore  is  desirable  to  develop  a  methodology  for  estimating  cladding 
panel  design  loads  from  strength  data  obtained  by  testing  small  standard 
specimens.     The  purpose  of  this  report  is  to  propose  such  a  methodology.  The 
methodology  utilizes  recent  advances  [3,  14] ,  which  enable  the  probability 
distribution  of  the  load  capacity  of  glass  panels  to  be  estimated  by  numerical 
methods  on  the  basis  of  information  on  the  strength,  stiffness,  and  fracture 
mechanics  properties  of  glass. 

Information  characterizing  the  elastic  behavior  of  glass  is  available  or  can 
be  obtained  by  well  established  test  methods.     The  elastic  parameters  exhibit 
relatively  small  variability,  so  that  uncertainties  pertaining  to  their  actual 
values  do  not  significantly  affect  the  estimation  of  design  loads.     In  recent 
years  test  methods  have  been  developed  to  obtain  fracture  mechanics  parameters 
governing  subcritical  crack  growth  and  strength  degradation  under  load.  Assoc- 
iated with  the  fracture  mechanics  parameters  are  large  uncertainties,  which 
include  uncertainties  with  respect  to  the  shape  of  surface  flaws.     However,  it 
is  shown  in  this  work  that  the  consequences  of  these  uncertainties  are  minor 
from  an  engineering  point  of  view. 

The  critical  question,  then,  is  whether  it  is  possible  to  obtain  from  tests  of 
small  standard  specimens  the  information  on  the  strength  of  glass  needed  to 
estimate  reliably  the  load  capacity  of  cladding  panels  whose  failure  is  due  to 
surface  flaws.     Evidence  presented  in  this  work  strongly  suggests  that,  if  the 
ring-on-ring  testing  method  is  employed,  the  answer  to  this  question  is 
affirmative. 
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2.     ESTIMATION  OF  LOAD  CAPACITY  OF  CLADDING  PANELS 


The  purpose  of  this  section  is  to  review  pertinent  developments  that  make  it 
possible  to  estimate  numerically  the  probability  distribution  of  the  load 
capacity  of  cladding  panels  if  the  fundamental  parameters  characterizing  the 
behavior  of  glass  are  known. 

2.1     FRACTURE  MECHANICS  OF  GLASS 


The  basic  criterion  for  fracture  is  derived  from  the  Griffith  equilibrium 
expression,  and  may  be  written  as 

Ki  =  Kic  (1) 


where  Kj  =  stress  intensity  factor,  and  Kjq  =  critical  value  of  Kj.  If 
equation  1  holds,  the  system  reaches  the  state  of  instability  wherein  the  rate 
of  crack  growth  becomes  for  practical  purposes  infinite  and  failure  occurs  [8, 
9].     Kj^Q  is  a  property  of  the  material  and  is  determined  experimentally.  The 
stress  intensity  factor,  Kj,  is  proportional  to  the  actual  stresses  in  the 
material  in  the  presence  of  cracks  causing  stress  concentrations.     Kj  can  be 
expressed  as  follows  [8,  9]: 

Ki(t)  =  Ya(t)  v^TcIF  (2) 


where 


Y  =  geometric  shape  factor, 

a  =  nominal  stress  (i.e.,  stress  calculated  by  assuming  the  absence  of 
cracks) , 

c  =  length  of  crack  normal  to  the  stress 
t  =  time. 


The  geometric  shape  factor,  Y,  in  equation  2  is  a  function  of  crack  geometry 
and  plate  dimensions  [10].     For  semi-elliptical  surface  flaws  (figure  1)  for 
which  c/h  =  0  and  a/b      0,  values  of  Y  calculated  in  reference  10  for  various 
ratios  c/a  are  shown  in  table  2.1. 

According  to  experiments  reported  in  references  11,  12,  and  13,  the  following 
relationship  holds  for  the  rate  of  subcritical  crack  growth  in  the  region 
Ki  <  Kic: 


lf=4K?(t)  (3) 

The  parameters*^  and  n  depend  upon  ambient  humidity  and  temperature  and  are 
obtained  experimentally.     Equation  3  expresses  quantitatively  the  fact  that  the 
cracks  in  an  element  of  glass  subjected  to  stress  for  some  length  of  time  will 
grow  -  albeit  not  catastrophically  -  provided  that  the  stress  is  contained 
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Table  2.1    Dependence  of  Factor  Y  Upon  Ratio  c/a  for  c/h      0  and  a/b  =  0 


c/ a 

Y 

0.00 

1.985 

0.25 

1.827 

0.50 

1.581 

0.75 

1.353 

1.00 

1.163 

1.10 

1.157 

1.20 

1.149 

1.30 

1.139 

1.40 

1.127 

1.50 

1.116 

1.60 

1.103 

1.70 

1.090 

1.80 

1.077 

1.90 

1.064 

2.00 

1.050 

within  a  certain  range.     This  phenomenon  is  referred  to  as  static  or  dynamic 
fatigue  according  to  whether  the  stress  is  constant  or  time-dependent. 

It  follows  from  equations  1  and  2  that  the  strength  of  glass,  S,  i.e.,  the 
value  of  the  nominal  stress  at  which  failure  occurs,  is 

^IC 

S(t)   (4) 

Y/^(tT 

If  Kj  and  c  are  eliminated  from  equations  2,  3,  and  4  and  the  notation  S(0)  = 
Sj^  is  used  (S-^  =  initial  strength),  the  following  relationship  is  obtained: 


S(t)  =  [S^  ^  -  1  /  a"  (t)  dx]      ^  (5) 


B  0 


where 


1/B  =  -^Ay2  (5a) 

[11].  If  we  consider  an  area,  A\^,  over  which  the  tension  stress,  a(t),  is  uniform 
and  independent  of  direction,  failure  occurs  within  that  area  if 

a(t)  >  S(Ak,  t)  (6) 
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where  S(A]^,  t)  =  strength  calculated  by  equation  5  in  which  Sj^  corresponds  to 
the  largest  initial  flaw  within  the  area  A^^,  regardless  of  the  direction  of 
that  flaw,  that  is 


(7) 


and  CjnaxC-^k,        ~  length  of  largest  flaw  within  Aj^  at  time  t  =  0. 


We  now  consider  the  case  in  which  the  state  of  stress  is  uniform  over  the  area 
A]^,  but  the  principal  stresses  within  the  area  are  unequal.     The  effect  of 
normal  stress  is  by  far  the  strongest  as  far  as  crack  propagation  is  concerned 
[9,  p.  54],  and  that  the  effect  of  shear  stresses  may  therefore  be  neglected 
(see,  e.g.,  reference  14).     The  center  of  A^  is  denoted  by  M]^,  and  the  normal 
tension  stresses  at  time  t  parallel  to  direction        are  denoted  by  aCM^,  ajj,,  t). 
We  define  the  sector  Aaj^,,  centered  on  the  direction  a^,  and  such  that  normal 
stresses  at  time  t  parallel  to  any  radial  direction  contained  within  that 
sector  may  be  assumed  to  differ  negligibly  from  a(Mk,  aj^,,  t)   [figure  2]. 
Failure  within  the  area  Ay^  will  not  necessarily  be  initiated  by  a  flaw  normal, 
or  almost  normal,  to  the  largest  principal  stress,  since  the  largest  of  these 
flaws  may  well  be  relatively  small.    Neither  will  failure  be  necessarily  ini- 
tiated by  the  largest  flaw  within  the  area  A^,  since  that  flaw  may  well  be 
perpendicular  to  a  relatively  low  normal  stress.    Rather,  failure  will  be  ini- 
tiated by  the  largest  of  the  flaws  normal  to  any  radius  within  the  sector  Aajj, 
[whose  length  is  denoted  by  c^-^  (A^,  Aa^,  t)]  to  which  there  corresponds  a 
strength,  S(A]^,  Aa^,  t),  such  that: 


For  any  sector,  Aa^,  the  strength  at  time  t,  S(A]^,  Aa^,  t),  is  calculated  by 
equation  5  in  which    ct(Mj^,  a^,  t)  is  substituted  for  aCx)  and  the  initial 
strength  S^CA^^,  Aa^)  is  substituted  for  Sj^.     The  following  relation  holds: 


where  Cu^x^^k*  Aajj,,  0)  =  length  of  largest  initial  crack  normal  to  any  radius 
within  the  sector  Aaj2,. 

2.2     RELATION  BETWEEN  INITIAL  STRENGTH  AND  THE  FINAL  STRENGTH  CORRESPONDING  TO 
A  60-SEC  LOAD 

A  particular  case  of  practical  interest  is  that  in  which  the  load  acting  on  the 
panel  has  constant  value  over  a  60-sec  time  interval  and  is  zero  outside  that 


t)  >  S(Ai^,  Aa£,  t) 


(8) 


S^(A^^,  Aaji)  = 


1/2 

^<^max(^k'  0) 


(9) 
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interval^.     The  60-sec  load  induces  at  point  Mj^  a  normal  stress  parallel  to 
the  direction  a^,  denoted  by  a5o(%>  ot^),  which  is  constant  throughout  the 
duration  of  the  load  and  equal  to  zero  at  all  other  times.     Failure  initiated 
by  a  flaw  normal  or  almost  normal  to  a  radius  within  the  sector  Aaji  occurs  at 
time  t  =  60  sec  if  the  stress  <760(^k>  equal  to  the  strength  SCAj^,  Aajj,, 

t  =  60).     Substituting  <75o(Mk»  <^0  <^it)  and  S^CAj^,  Aaj^,)  for        in  equation  5, 

it  follows  that  the  failure  stress  is  given  by  the  relation 


^60 


(Mk,  a^)  [olo  (Ml,,  a^)  +  |q ]  +  |q  sJ  ^Aj,,  Aa;^)  (10) 

B/MD^  ^2 


For  soda-lime  glass  ^^^Q^^y^i  ot^)  is  of  the  order  of  10"^(MPa)  or  more,  and  B/60 
is  of  the  order  of  l(MPa)     or  less,  so  that  the  failure  stress  may  be 


written  as 

n-2 


^i^  (A^,  Aaj^) 


%0(%»  ajl)  =^   TT —  (11) 

1  1/n 

(60  i) 

2.3     ESTIMATION  OF  60-SEC  LOAD  CAUSING  FAILURE  OF  A  PANEL 

The  relationship  between  the  load  acting  on  a  cladding  panel  and  the  stresses 
in  the  panel  is  generally  nonlinear.     This  relationship  can  be  obtained  by 
using,  for  example,  a  finite-difference  program  such  as  that  developed  by 
Texas  Tech  University  [15],  and  is  commonly  expressed  in  terms  of  the  nondimen- 
sional  quantities 

LF  =  pb^/Dh  (12) 

SF  =  ab^h/D  (13) 

where  LF  and  SF  =  loading  and  stress  factors,  respectively,  p  =  uniform  load 
per  unit  area,  a  =  stress,  b  =  smaller  side  of  rectangular  plate,  D  =  flexural 
rigidity,  defined  by: 

D  =  (14) 
12(l-v2) 

where  E  =  modulus  of  elasticity,  v  =  Poisson's  ratio,  and  h  =  thickness  of  the 
plate. 


^  The  60-sec  load  is  the  standard  reference  load  used  in  glass  cladding  design  charts 
in  the  United  States. 
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Once  the  relationship  between  loads  and  stresses  is  known,  it  is  possible  to 
obtain  the  60-sec  loads,  P^q*  corresponding  to  the  failure  stresses  o^q  (Mj^,  a^) 
calculated  by  equation  11  for  each  point  }\  of  a  sufficiently  dense  grid  and 
for  each  of  a  sufficient  number  of  directions  aj^.     The  smallest  of  these  loads, 
denoted  by  P50>  is  the  load  causing  failure  (i.e.,  the  load  capacity)  of  the 
panel  characterized  by  the  set  of  initial  strengths  Sj[(A]^,  Aajj).  Calculations 
of  the  load  capacity,  P60>  carried  out  for  a  large  number,  M,  of  panels. 

The  probability  distribution  of  60-sec  load  capacity  can  be  estimated  from  the 
M  values  so  obtained.     A  computer  program  for  estimating  this  distribution  if 
the  initial  strengths  S-|[(A;^,  Aajj)  are  specified  is  described  in  Appendix  III. 
In  the  case  of  heat-strengthened  or  tempered  glass,  the  procedure  for  estimating 
the  probability  distribution  of  the  60-second  load  capacity  p^g  is  the  same, 
except  that  stresses  cy5o(Mic>0'e^~'^R  should  be  used  in  equation  11  in  lieu  of  the 
stress,  0()QiKy^,ai)  ^  where  Or  denotes  the  residual  thermal  stress,  which  can 
be  determined  by  routine  experimental  procedures. 

2.4     SPECIFICATION  OF  INITIAL  STRENGTHS  SjCA]^,  Aaj^) 

The  initial  strength  Sj^(A^)  corresponding  to  the  largest  flaw  within  Aj^, 
regardless  of  its  direction,  is  commonly  described  by  a  Weibull  distribution: 

Si(Ak)  -  Ps  ™ 

P[Si(Ak)]  =  1  -  exp  {  -[       g^(^^)     ]   }  (15) 

From  the  assumption  that  the  number  of  flaws  of  any  given  size  is  on  the  average 
proportional  to  the  area  A  being  considered,  it  follows  that 

Si(A)  -  Pg  m 

P[Si(A)]  =  1  -  exp  {  -[  — s-(^)  ]   }  ~  (16) 


where 


m 


S^(A) 


=  (- 


(16a) 


[25,  pp.  5  and  10].  We  will  refer  to  equation  16  as  the  fundamental  Weibull 
distribution  of  the  strength  of  glass. 
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Similarly,  from  the  assumption  that  the  flaw  orientations  are  uniformly 
distributed  [i.e.,  that  the  number  of  flaws  normal  to  the  stress  o^o^^k'  ^l^^ 
is  on  the  average  equal  to  Aajj,/(iT/2 )  times  the  number  of  flaws  parallel  to  any 
direction],  it  follows  that  the  probability  distribution  of  the  initial  strength 
Si(Ak,  Aaji)  is 


S^CAj^    Aaj^  )  ~  l^s 

P[S,(A^.  L.,n  =  1  -  exp  1  -l^J--^^  r]  (17) 

where 


1/m 

S^CA^,  Aaji)  =  (2^)        S^CAk)  (18) 

Assuming  that  equations  15  through  18  are  valid,  it  follows  that  the  initial 
strength  S-j^(A]^,  Aajj)  can  be  specified  by  the  probability  distribution  given  by 
equation  17,  provided  that  the  parameters  Sq(A),   yg,  and  m  are  known. 
These  parameters  can  be  obtained  from  ring-on-ring  tests,  which  are  described 
in  the  following  section. 


^  Or  to  the  stress  o^qQA^,  -aj^)  =  ct^qCMj^,  a^^). 
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3.     RING-ON-RING  TESTS 

3.1     PRINCIPLE  AND  DESCRIPTION  OF  RING-ON-RING  TESTING 


Ring-on-ring  testing  devices  involve  the  creation  of  a  state  of  uniform 
axi symmetric  tension  stress  in  the  central  portion  of  one  of  the  faces  of  a 
circular  plate.     This  can  be  accomplished  by  placing  the  plate  on  a  segmented 
circular  ring  and  by  applying  on  its  upper  surface  a  load  transmitted  through  a 
circular  ring  concentric  with  and  having  a  smaller  diameter  than  the  segmented 
support . 

As  noted  in  reference  17 ,  because  of  the  high  elastic  modulus  and  hardness  of 
glass,  any  nonperfect  contact  between  rigid  loading  rings  and  the  test  sample 
can  lead  to  deviations  from  axisymmetry  in  the  stress  field.     Reference  17 
describes  a  ring-on-ring  device  designed  to  eliminate  such  nonunif ormities. 
Each  ring  consists  of  a  closely  wound  coil.     The  load  that  each  coil  transmits 
to  the  plate  is  applied  by  a  rubber  diaphragm  which  covers  a  circular  groove 
filled  with  fluid.     The  function  of  the  fluid  is  to  equalize  the  loading  along 
the  coils.     Strain  gage  measurements,  and  measurements  of  strength  of  indented 
specimens,  are  reported  in  reference  17  for  ringon-ring  devices  both  of  the 
rigid  type  and  of  the  type  just  described.     According  to  the  results  of  refer- 
ence 17,  errors  in  the  measurement  of  strengths  due  to  the  use  of  rigid  ring- 
on-ring  devices  are  less  than  5  percent  for  annealed  soda-lime  float  glass. 
The  errors  are  considerably  larger  (about  20  percent)  for  thermally  tempered 
crown  glass. 

The  uniform  stress  in  the  central  portion  of  the  tension  face  of  a  circular 
plate  subjected  to  a  ring-on-ring  loading  test  is 

a  =  -1^  [2(1  +  V)  An  ^  +    (l-v)(a2-b2)  ajj 
4irh2  b  a2  r2 

where  p  =  load,  h  =  plate  thickness,  a  =  radius  of  the  support  ring,  b  =  radius 
of  the  loading  ring,  R  =  radius  of  the  disc,  and  v  =  Poisson's  ratio  [18]. 

For  the  case  b/a  =  1/2,  the  distribution  of  the  stress,  a,  along  a  radius,  r, 
is  represented  in  figure  3  [19],     Note  that,  owing  to  local  effects,  the 
central  portion  of  the  plate  within  which  the  stress  a  is  given  by  Eq.  19  is 
defined  by  the  relation  r  ^  0.705b, 

It  follows  from  appendix  I  that  the  precision  of  strength  estimates  based  on  the 
testing  of  any  given  sample  of  ring-on-ring  specimens  Increases  with  the  number 
of  specimens  within  that  sample  that  have  failure  origins  in  the  central, 
uniform  stress  area  of  the  specimen.     It  would  therefore  be  desirable  to  attempt 
the  development  of  test  devices  aimed  at  inhibiting  stress  corrosion  on  the 
tension  side  of  the  nonuniform  stress  area  of  the  specimen.     This  could  be 
accomplished  by  immersing  that  side  in,  or  covering  it  with,  an  inert  agent 
such  as  water-free  oil. 
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3.2     ESTIMTION  OF  INITIAL  STRENGTHS  FROM  RING-QN-RING  TEST  RESULTS,  AND 


INFLUENCE  OF  UNCERTAINTIES  WITH  RESPECT  TO  PARAMETERS  B  AND  n  UPON  THE 


ESTIMATION  OF  DESIGN  LOADS 

Ring-on-ring  tests  yield  the  stress  at  the  time  of  failure,  a(tf),  which  is 
equal  to  the  strength  at  the  time  of  failure,  S(tf).     To  estimate  the  initial 
strength  corresponding  to  S(tf),  equation  5  is  used.     In  the  case  of  a  ramp 
loading  on  a  specimen  that  exhibits  a  linear  load-stress  relationship  up  to 
failure  [21,  22],  steps  similar  to  those  that  led  to  equation  11  yield  the 
relation 

1/n  n-2 

a(tf)  -  [— ^iip-]        Si  (20) 
''f  B 

Equation  20  allows  the  estimation  of  Sj^  from  measurements  of  a(tf)  and  tf.  The 
estimates  of  S^  depend  upon  the  values  assumed  for  n  and  B. 

We  now  consider  a  specimen  subjected  to  a  ramp-like  load,  for  which  the  time  to 
failure  and  the  stress  at  the  time  of  failure  were  found  to  be  tf  and  a(tf), 
respectively.    We  seek  the  stress,  059 »  induced  by  a  constant  load  with  a  60-sec 
duration  that  would  have  caused  failure  of  that  same  specimen.     It  follows  from 
equation  11  (with  Aa^  =    it/2)  and  equation  20  that 

^60==  — r;  ^^^f^  ^21) 

[(n+l)  60 
tf 

Equation  21  shows  that  the  ratio  a(^Q/aitf)  is  independent  of  B.     Thus,  even 
though  the  variability  of  B  is  fairly  large,  reflecting  as  it  does  uncertainties 
with  respect  to  ^jq,  n.  A,  and  Y  (see  equation  5a  and  table  2.1),  the  effect  of 
this  variability  upon  the  ratio  0(^Q/aitf)  can  be  ignored. 

Table  3.2  lists  results  of  calculations  which  show  that  for  ratios  tf/60  of  the 
order  of  0.5  to  2  the  effect  of  the  parameter  n  upon  the  ratio  o^^q/ o(tf) 
is  small^.     The  calculations  were  based  on  two  values:     n=16  and  n=19.7 
[23].     These  values  are  consistent  with  results  of  tests  conducted  within  the 
framework  of  this  project  on  indented  soda-lime  specimens,  in  accordance  with 
the  method  described  in  reference  24.     The  tests  yielded  the  value  n=15.95 
when  measurement  results  were  not  corrected  for  residual  stress  effects  at  the 
indented  crack  tips,  and  the  value  n=20.6  when  the  residual  stress  effects 
were  accounted  for.     The  magnitude  of  the  residual  stresses  at  the  crack  tips 
of  a  glass  panel,  and  therefore  the  effective  value  of  n  that  should  be  used 
in  determining  fatigue  effects  for  the  panel,  is  unknown.     However,  it  is 
reasonable  to  assume  that  the  latter  will  lie  between  the  values  n=15.95  or  so 


^  The  choice  of  a  proper  loading  rate  can  in  practice  ensure  that  the  ratios  t^/60 
are  indeed  of  that  order  -  see  tables  3.3  and  3.4. 
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and  n=20.6  or  so.  Hence  the  choice  of  the  values  n=16  and  n=19.7  for  the 
sensitivity  calculations  of  table  3.2. 

Table  3.2    Ratios  a^^Q/oitf)  for  tf  =  30  sec  and  tf  =  120  sec, 
Assuming  n  =  19.7  and  n  =  16.00 


n  =  19.7 

n  =  16 

Difference 

tf  =  120  sec 

0.888 

0.875 

1.5% 

tf  =  30  sec 

0.828 

0.802 

3% 

Equation  21  and  table  3.2  (which  show,  respectively,  that  the  ratio  o^^q/ oitf) 
Is  Independent  of  B,  and  that  for  any  given  stress,  a(tf),  obtained  by  testing 
a  rlng-on-rlng  specimen,  the  corresponding  60-sec  strength,  059,  depends  weakly 
upon  n)  suggest  that  estimates  of  the  probability  distribution  P(p5q),  of 
the  60-sec  load  capacity  of  a  panel,  P6o»  Inferred  on  the  basis  of  results 
of  rlng-on-rlng  tests  are  also  Independent  of  B  and  weakly  dependent  upon  n. 
This  was  confirmed  by  estimates  of  PCp^q)  carried  out  for  a  4  ft  x  4  ft  x 
1/8  in.  annealed  glass  panel  in  which  several  assumed  sets  of  values  of  n  and 
1/B  were  used  to  convert  strengths  at  time  of  failure,  a(tf),  into  initial 
strengths, 

3.3     RESULTS  OF  RING-ON-RING  TESTS 

This  section  presents  results  of  ring-on-ring  tests  conducted  within  the 
framework  of  this  project.     The  ring-on-ring  testing  device  (figure  4)  consisted 
of  rigid  rings  with  radii  a  =  0.0603  m  (support  ring)  and  b  =  0.0254  m  (loading 
ring).     The  device  was  employed  in  combination  with  a  lOK  ibf  Universal  Testing 
Machine.    All  the  specimens  were  subjected  to  ramp  loads  (i.e.,  loads  Increasing 
for  practical  purposes  linearly  with  time).     The  glass  used  in  the  tests  was  new 
and  was  obtained  from  the  same  manufacturer  and  batch. 

Tests  were  conducted  in  air  on  a  set  of  56  annealed  float  glass  square  specimens 
with  side  D  =  7  in.   (0.1792  m)  and  on  a  set  of  29  annealed  float  glass  circular 
specimens  with  radius  R  =  3.5  in  (0.889  m) .     The  nominal  thickness  was  h  =  1/4 
in.  (6  mm)  for  all  specimens.     The  stresses  at  the  center  of  the  plates  were 
calculated  by  equation  19,  in  which  h  was  the  measured  (rather  than  the  nominal) 
thickness  for  each  specimen.     Following  reference  19,  in  the  case  of  the  square 
plates  the  parameter  R  in  equation  19  was  assumed  to  be  equal  to  one-half  the 
average  of  the  edge  and  diagonal  lengths.     This  was  confirmed  experimentally 
to  within  a  few  percentage  points  by  strain-gage  measurements  of  stresses  on 
circular  and  square  rlng-on-rlng  specimens.     The  test  results  and  the  calculated 
values  of  the  stresses  at  the  center  of  the  plates  at  the  time  of  failure,  0(tf), 
are  listed  in  tables  3.3  and  3.4.     It  is  noted  that  the  influence  of  humidity 
upon  the  results  of  tables  3.3  and  3.4  is  negligible  (see  reference  14,  p. 
35). 
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If  it  is  assumed  that  n  =  19.7  and  1/B  =  0.0738569(MPa)~2s~l  (corresponding  to 

j»  n 

<A0  =  1.08  (MPa)""        2  s"^  [23],  Kj^,  =  0.75  MPa  [23],  and  Y  =  1.12  -  see  equa- 
tion 5a),  it  follows  from  equation  20  that  the  initial  strengths  are 

=  0.7272  t^l^''-'^  [a(tf)]l^-^/l^-^  (22) 

Statistical  analyses  of  the  initial  strengths  calculated  by  equation  22  from 
the  stresses  a(tf)  were  carried  out  as  indicated  in  appendix  I.     It  was 
assumed  that  the  area  A  referred  to  in  appendix  I  is  a  circle  with  radius 
r  ^  0.705b  =  0.705  in.  (0.0179  m) ,  i.e.,  A  =  TTr2      1.56  in. 2  (0.001  m2) 
(see  figure  3). 

Estimated  parameters  of  the  Weibull  distributions  of  the  initial  strengths,  Sj^, 
from  results  of  tests  used  in_conj unction  with  equation  22  are  listed  in  table 
3.5,  where  the  sample  means,  S-j^,  standard  deviations,  s(Sj^),  coefficients  of 
variation,  s(S£)/Sj^,  sample  maximum,  Sj_jjjg2j,  and  sample  minimum,  Sj^^j^j^, 
of  the  data  are  also  shown. 

The  initial  strengths  Sj^  calculated  by  equation  22  are  nominal,  rather  than 
actual,  since  the  values  of  B  and  n  used  therein  are  uncertain.     However,  it 
follows  from  the  form  of  equations  11  and  20  and  from  the  results  of  table  3.2 
that,  if  the  values  of  B  and  n  used  in  these  equations  are  the  same,  the  effect 
of  uncertainties  with  respect  to  the  actual  values  of  B  and  n  largely  cancels 
out  when  estimating  the  load  capacity  of  glass  panels. 
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Table  3.3    Test  Data  for  Circular  Plates 
h  tf  p  r  RH%  a(tf) 


1 

0.5334 

45 

2358.5 

0.76 

67 

47.32 

2 

0.5639 

54 

3453.2 

0.76 

67 

61.99 

3 

0.5690 

59 

3419.9 

1.27 

71 

62.06 

4 

0.5613 

48 

2839.1 

1 .78 

71 

51.44 

5 

0.5537 

73 

4534.6 

1.27 

71 

84.44 

6 

0.5613 

43.5 

2109.3 

2.03 

66 

38.22 

7 

0.5080 

36 

1637 .6 

1.27 

66 

36.23 

8 

0.5080 

41.5 

2073.7 

2.29 

66 

45.87 

9 

0.5334 

31 

1713.3 

2.03 

60 

34.38 

10 

0,5486 

82 

4694.8 

2.29 

60 

89.05 

11 

0.5436 

57 

2536.5 

0.00 

67 

49.00 

12 

0.5512 

59 

3924.9 

1 .27 

67 

73.75 

13 

0.5334 

47 

3021 .6 

1 .02 

67 

60.63 

14 

0.5359 

70 

3742.5 

1 .52 

67 

74.39 

15 

0.5537 

45 

2460.9 

1.27 

72 

45.82 

16 

0.5588 

52 

2963.7 

1.27 

72 

54.18 

17 

0.5563 

45 

2336.3 

2.29 

72 

43.10 

18 

0.5563 

73 

3898.2 

2.03 

70 

71.91 

19 

0.5537 

51 

2269.5 

1.52 

70 

42.26 

20 

0.5537 

50 

2216.1 

2.29 

70 

41.26 

21 

0.5461 

73 

3804.8 

1.52 

70 

72.83 

22 

0.5486 

54 

2594.4 

2.03 

70 

49.21 

23 

0.5359 

71 

3675.7 

0.51 

70 

73.07 

24 

0.5385 

76 

5255.5 

2.54 

60 

103.46 

25 

0.5410 

76 

3960.5 

2.54 

72 

77.25 

26 

0.5537 

62 

3956.1 

2.79 

71 

73.66 

27 

0.5461 

74 

4000.6 

2.79 

67 

76.58 

28 

0.5639 

49 

2523.2 

3.30 

66 

45.30 

29 

0.5334 

37 

2055.9 

3.30 

67 

41.25 

where         h  =  thickness  in  cm 

tf  =  load  time  in  seconds 
p  =  failure  load  in  Newtons 

r  =  distance  from  origin  of  fracture  to  center  in  cm 
RH  =  relative  humidity 
oXtf)  =  final  strength  in  MPa 


12 


Table  3.4     Test  Data  for  Square  Plates 


h 

f-  £ 

r 

1 

0.5385 

2647 . 8 

D  / 

9  9Q 

7  /< 

en  no 

2 

0.5359 

2287  .3 

54 

n  7  fx 

IL 

'♦J  .03 

3 

0.5563 

2429.7 

48 

0  7f» 

U  •  /  D 

1 L 

AT  n7 

4 

0.5537 

2621 . 1 

61 

2 .03 

\j\> 

AA  Qfl 

HO  .  7U 

5 

0.5461 

2242.8 

49 

1  7Q 

DO 

1  .  ZO 

6 

0.5461 

391 1 .6 

73 

79 

7  1  OS 

7 

0.5461 

3773 .6 

81 

79 

AQ   A  9 

07  mHL 

8 

0.5436 

2883.6 

74 

1  97 

79 

9 

0.5461 

4427.8 

95 

0.76 

66 

9<\  AS 

10 

0.5461 

3017 . 1 

73 

0.51 

66 

11 

0.5410 

2451 .9 

67 

1  7Q 

A9 
DZ 

AS  QA 
J  •  70 

12 

0.5334 

4561  .3 

100 

7n 

R7   Q  S 

13 

0.5334 

4525 . 7 

1  nn 

1  70 

/  u 

ft7  9(i 
o  /  .  ZD 

14 

0.5334 

3163.9 

80 

1  97 

70 

O  i  •  uu 

15 

0.5436 

2309.6 

63 

1 .02 

A7  8R 

16 

0.5359 

2919.2 

75 

2 .29 

62 

SS  76 

17 

0.5537 

2558.8 

67 

1.27 

66 

AS  7Q 

18 

0.5461 

5157  .6 

97 

1 .27 

66 

94 .87 

19 

0.5486 

5736. 1 

106 

1.52 

66 

1  OA  SA 

20 

0.5436 

5006.3 

101 

2.29 

62 

92.94 

21 

0.5410 

3840.4 

73 

1.52 

62 

71 .98 

22 

0.5410 

4583.5 

93 

2.29 

62 

85.91 

23 

0.5349 

3119.5 

73 

1.52 

62 

59.59 

24 

0.5334 

4401 .1 

93 

2 .03 

68 

84  8S 

25 

0.5512 

3408.7 

84 

2.03 

68 

61 .55 

26 

0.5410 

5402 .3 

104 

2 . 29 

62 

1  n  1  96 

27 

0.5512 

4458.9 

93 

0.76 

62 

80.51 

28 

0.5359 

4579.1 

103 

1 . 27 

62 

87  .47 

29 

0.5588 

4570. 2 

92 

1.02 

62 

80.29 

30 

0.5512 

1984.7 

51 

1.79 

66 

35.84 

31 

0.5461 

4423.3 

98 

2.03 

66 

81.37 

32 

0.5436 

5816.2 

117 

2.54 

66 

101 .98 

33 

0.5385 

2558.8 

56 

2.54 

62 

48.41 

34 

0.5512 

3399.8 

82 

2.54 

70 

61.39 

35 

0.5410 

4783.8 

103 

2.54 

62 

89.66 

36 

0.5512 

2656.7 

73 

2.54 

68 

47 .97 

37 

0.5461 

5388.9 

102 

2.54 

68 

99.13 

38 

0.5537 

5647 . 1 

107 

2.54 

68 

101 .05 

39 

0. 5461 

5660.4 

111 

2.54 

68 

104. 12 

40 

0.5436 

4147 .4 

89 

2.54 

66 

76.99 

41 

0.5563 

5024 . 1 

101 

2.54 

62 

89.06 

42 

0.5385 

4104.2 

101 

2.54 

62 

96.56 

43 

0.5334 

2233.9 

65 

2.54 

62 

43.07 

44 

0.5385 

5526.9 

110 

2.54 

68 

104.56 

45 

0.5410 

2710.1 

55 

2.79 

66 

50.80 

46 

0.5410 

3088 .3 

63 

2.79 

66 

57.89 

47 

0.5359 

3964.9 

94 

2.79 

62 

74.74 

48 

0.5385 

3333.1 

84 

2.79 

70 

63.06 

49 

0,5461 

5611.5 

114 

2.79 

70 

103.22 

50 

0.5563 

5117.5 

100 

2.79 

62 

90.72 

51 

0.5359 

4699.2 

87 

3.05 

66 

89.76 

52 

0.5358 

4294.3 

97 

3.05 

62 

82.03 

53 

0.5486 

3288.6 

80 

3.05 

62 

59.94 

54 

0.5334 

2621.1 

74 

3.30 

68 

50.54 

55 

0.5486 

2843.6 

74 

3.56 

68 

51.83 

56 

0.5334 

4325.4 

97 

3.56 

62 

83.40 

where 

h  =  thickness  in  cm 

tf  =  load  time  in  seconds 
p  =  failure  load  in  Newtons 

r  =  distance  from  origin  of  fracture  to  center  in 
RH  =  relative  humidity 
a(tf)  =  final  strength  in  MPa 
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4.     ESTIMATES  OF  THE  LOAD  CAPACITY  OF  GLASS  PANELS 


The  methodology  for  estimating  the  60-sec  load  capacity  of  glass  panels 
proposed  in  this  work  consists  of  using  the  procedure  described  in  section  2.3 
in  conjunction  with  equations  17  and  18.     The  parameters  SqCA]^),  m,  and  Pg  in 
these  equations  are  obtained  by  fitting  a  Weibull  distribution  to  the  nominal 
initial  strengths  Sj^  calculated  from  the  breaking  strengths  a(tf)  by  using 
equation  20  (or  a  similar  relation  if  the  dependence  of  load  on  time  is  not 
linear).     The  values  of  the  parameters  B  and  n  used  in  equation  20  must  be 
the  same  as  those  of  equation  11.     Values  of  n  used  in  the  calculations  should 
be  based  on  experimental  results  obtained,  e.g.,  by  techniques  described  in 
reference  24. 

This  methodology  was  applied  to  a  4  ft  x  4  ft  x  1/8  in  (1.22  m  x  1.22  m  x  3  ram) 
annealed  glass  panel  simply  supported  on  four  sides,  using  the  values  E  = 
68.9  GPa,  V  =  0.22,  the  fracture  mechanics  parameters  listed  in  section  3.3, 
and  the  parameters  of  the  Weibull  distribution  of  the  initial  glass  strength 
listed  in  table  3.5.     The  grid  size  and  the  angle  Aa£  used  in  the  numerical 
calculations  were  7.62  mm  x  7.62  mm  and  18°,  respectively.     For  each  of  the 
sets  of  two  and  three  parameters  listed  in  table  3.5  the  load  capacities,  P60> 
of  1000  of  panels  were  estimated,  and  the  values  so  obtained  were  fitted  by 
two-  and  three-parameter  Weibull  distributions,  respectively.     The  parameters 
of  the  best  fitting  distributions  are  listed  in  table  4.1,  which  also  lists 
mean  values,  P60 >  standard  deviations,  s(p5o),  coefficients  of  variation, 
s(p5o)/p60»  lo^ds  corresponding  to  a  probability  of  failure  of  8  in  1000, 
P50 (0*008),  and  loads  corresponding  to  a  probability  of  failure  of  0.5. 

Estimates  of  the  load  capacity  corresponding  to  a  probability  of  failure  of  8  in 
1,000,  p5o(0.008),  based  on  full-scale  measurements  are  provided  for  design 
purposes  in  references  1  and  3.     According  to  reference  1,  p5o(0.008)  =  26  psf 
(1  psf  =  47.9  Pa).     According  to  reference  3,  P5o(0.008)  -  23  psf.     It  is 
noted  that  to  account  for  strength  degradation  in  service,  the  value  of  refer- 
ence 3  corresponds  to  strengths  reduced  by  a  factor  of  2/3  with  respect  to 
those  obtained  in  new  glass  [26].     Had  this  reduction  not  been  effected,  i.e., 
had  the  strength  of  new  glass  been  used,  to  the  value  P5o (0.008)  -  23  psf 
there  would  have  corresponded  roughly  the  value  P5q( 0.008)  =  (3/2)  x  23  = 
34.5  psf. 

The  values  p^gCO.OOS)  estimated  in  this  report  on  the  basis  of  three-parameter 
Weibull  distributions  (col.  9  of  table  4.1,  Case  I)  are  somewhat  higher  than 
the  corresponding  value  based  on  references  1  and  3  (35.5  psf  to  42.1  psf  versus 
26  psf  to  34.5  psf.)     These  differences  could  be  explained  by  two  factors. 
First,  the  estimates  obtained  in  this  report  do  not  take  into  account  edge 
failures,  which  can  reduce  the  load  capacity  of  panels  considerably.  (For 
example,  according  to  reference  7,  38  percent  of  the  total  number  of  panel 
failures  reported  therein  originated  at  the  edges.)     Second,  sampling  errors 
in  the  estimation  of  the  Weibull  parameters  of  the  strength  may  be  present, 
since  the  estimated  values  of  p5o(0.008)  are  based  on  the  testing  of  a  limited 
number  of  ring-on-ring  test  specimens  (from  29  for  case  1  to  85  for  case  3). 
Therefore,  as  shown  in  appendix  II,  the  true  values  of  p5o(0.008)  could  well 
be  as  low  as  about  20  psf  and  as  high  as  about  50  psf.     Note  also  that  the 
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estimated  coefficients  of  variation  of  the  load  capacity  p^Q  obtained  by 
using  three-parameter  Weibull  distributions  (0.09  to  0,105)  are  lower  than  the 
value  indicated  in  reference  3  (-  0.22).     On  the  other  hand,  as  can  be  seen 
from  table  4.2  which  summarizes  test  results  reported  in  reference  7,  coeffi- 
cients of  variation  of  the  load  capacity  of  new  glass  panels  obtained  from  any 
one  manufacturer  are  close  in  at  least  four  out  of  eight  cases  to  those  estimated 
herein. 

It  is  concluded  that  the  procedure  proposed  in  this  report  provides  credible 
estimates  of  the  load  capacity  of  panels  experiencing  surface  failures  provided 
that  three-parameter  Weibull  distributional  models  are  used.     The  precision  of 
these  estimates  depends  upon  number  of  ring-on-ring  test  results  on  which  they 
are  based,  as  suggested  by  numerical  experiments  presented  in  appendix  II. 

A  second  conclusion  of  interest  is  that  modeling  the  strength  of  glass  and  the 
load  capacity  of  cladding  panels  by  two-parameter  Weibull  distributions  leads 
to  a  clear  incompatibility  between  results  obtained  by  testing  ring-on-ring 
specimens  on  the  one  hand,  and  the  documented  behavior  of  cladding  panels  on 
the  other.     Indeed,  the  use  of  two-parameter  distributions  yields  estimates  of 
P^q(0.008)  of  the  order  of  1  psf  (see  table  4.1).     These  are  grossly  incompatible 
with  the  values  quoted  previously  from  references  1  and  3.    A  similar  conclusion 
was  reached  independently  by  Walker  [27]. 
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Table  4.1     Estimated  Statistics  and  Weibull  Distribution  Parameters  of  the 
1-min  Load  Capacity,  p^q  >  ^o^^  a  4  ft  x  4  ft  x  1/8  in  Annealed 
Glass  Panel  Supported  on  Four  Sides  (Based  on  Tests  in  Air) 


Case 

P60 
(psf) 

s(p6o) 
(psf) 

s(P60)/P60 

(P60)o 
(psf) 

1^60 
(psf) 

™&0 

P60(0.008) 
(psf) 

P60(0-5) 
(psf) 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

1 

42.8 

3.86 

0.090 

9.80 

34.14 

2.41 

35.5 

42.6 

2 

50.9 

5.00 

0.098 

11.04 

41.08 

2.05 

42.1 

50.3 

3 

45.4 

4.78 

0.105 

12.20 

34.57 

2.42 

36.2 

45.1 

1 

7.3 

3.43 

0.471 

8.23 

0 

2.26 

0.98 

7.0 

lib 

2 

9.5 

4.87 

0.515 

10.71 

0 

2.06 

1.03 

9.0 

3 

6.9 

3.58 

0.516 

7.86 

0 

2.06 

0.75 

6.6 

^  Based  on  three-parameter  Weibull  distributions  of  load  capacity  p^Q  and 
three-parameter  Weibull  distributions  of  the  strength  of  glass. 


Based  on  two-parameter  Weibull  distributions  of  load  capacity  p^Q  and 
two-parameter  Weibull  distributions  of  the  strength  of  glass. 

P50  ~  sample  mean 

s(P60^  ~  sample  standard  deviation 
(P50)o  ~  scale  parameter 
y^Q  =  location  parameter 
mgQ  =  shape  parameter 

P50(0«008)  =  1-min  load  capacity  corresponding  to  probability  of  failure  of 
panel  of  8  in  1,000 

P5q(0.5)  =  1-min  load  capacity  corresponding  to  probability  of  failure  of  0.5. 
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5.     SUMMARY  AND  CONCLUSIONS 


A  methodology  was  proposed  for  estimating  the  probability  distribution  of  the 
load  capacity  of  annealed  glass  panels  whose  failure  is  due  to  surface  flaws. 
The  methodology  requires  the  calculation  of  stresses  induced  on  the  panel 
surface  by  the  external  loads  and  employs  information  on  the  modulus  of  elas- 
ticity, Poisson's  ratio,  the  fracture  mechanics  parameters,  and  the  probabilis- 
tic description  of  the  glass  strength.     This  description  is  obtained  by  using 
the  ring-on-ring  testing  method.     It  is  shown  that  owing  to  the  way  in  which 
this  method  is  used,  errors  in  the  estimation  of  the  load  capacity  due  to 
uncertainties  with  respect  to  the  fracture  mechanics  parameters  of  glass  and 
to  the  shape  of  surface  flaws  are  largely  cancelled.     Results  of  ring-on-ring 
tests  and  of  calculations  based  on  those  tests  show  that  credible  predictions 
of  panel  load  capacities  are  obtained  that  are  consistent  with  data  available 
in  the  literature,  provided  that  the  probability  distribution  of  the  glass 
strength  and  the  probability  distribution  of  the  panel  load  capacity  are  modeled 
by  the  three-parameter  Weibull  distribution.     The  two-parameter  Weibull  distri- 
bution appears  to  provide  an  incorrect  model  of  the  strength  of  glass  and  of 
the  load  capacity  of  glass  panels. 

The  results  obtained  in  this  work  strongly  suggest  that  the  proposed  methodology 
can  provide  reliable  estimates  of  the  load  capacity  of  glass  panels  whose 
failure  is  due  to  surface  flaws.     However,  a  definitive  statement  to  this  effect 
would  require  validation  based  on  (1)  larger  ring-on-ring  test  samples  than 
those  used  in  this  work,  and  (2)  reliable  statistics  of  the  load  capacity  of 
panels  manufactured  from  the  same  batch  of  glass  as  the  ring-on-ring  test 
specimens . 

The  topic  of  glass  panel  failures  due  to  edge  imperfections  was  not  addressed 
in  this  report.     As  shown  by  the  statistics  in  table  4.2,  the  ratio  of  such 
failures  to  the  total  number  of  failures  is  highly  variable.     It  is  suggested 
that  the  probability  distribution  of  the  load  capacity  of  glass  panels,  regard- 
less of  type  of  failure,  can  be  modeled  from  information  on  statistics  of  edge 
failures  on  the  one  hand  and  of  surface  failures  on  the  other. 

In  view  of  their  relatively  low  cost,  ring-on-ring  tests  could  be  more  economical 
than  full-size  panel  tests,  in  spite  of  the  relatively  large  number  of  specimens 
that  would  have  to  be  tested  in  order  to  attain  acceptable  precisions  of  the 
estimates.     The  extent  to  which  this  is  the  case  would  have  to  be  determined 
by  studies  based  on  more  extensive  test  data  than  have  been  obtained  within 
the  framework  of  this  project. 

Ring-on-ring  testing  may  be  a  desirable  alternative  to  full-size  panel  tests 
not  only  for  economical  reasons,  but  also  in  situations  where  the  amount  of 
material  available  for  testing  is  limited.     This  might  be  the  case  in  studies 
of  in-service  strength  degradation  in  which  the  material  being  tested  consists 
of  weathered  window  glass  recovered  from  existing  buildings. 
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Figure  1.  Notations 


Figure  3.  Tangential  and  Radial  Stress  Distribution  (a 
b/a  =  1/2;  b  =  radius  of  loading  ring) 


radius  of  support  ring; 


APPENDIX  I.     MAXIMUM  LIKELIHOOD  ESTIMATION  FOR  THE  THREE  PARAMETER  WEI  BULL 
DISTRIBUTION  BASED  ON  A  SAMPLE  WITH  CENSORING  DUE  TO  COMPETING 
RISKS 


Let  X  be  a  random  variable  representing  the  tensile  strength  of  glass  cladding 
specimens  within  an  area  subjected  to  uniform  axial  bending.     It  is  assumed 
that  the  distribution  of  X  is  Weibull  that  is, 


where  a,  3>0,  X>0,  x>X  ,  and  a,  3,  and  X  are  the  three  parameters  of  the 
distribution. 

We  assume  that  a  number  of  specimens,  say  n,  are  tested  and  each  of  the  specimens 
has  the  same  Weibull  distribution  of  tensile  strengths  X.     However,  what  we 
observe  at  each  trial  is  either  the  tensile  strength  at  which  the  specimen 
fails  or  the  tensile  strength  outside  the  area  A.     For  a  gien  specimen,  let  Q 
denote  the  strength  at  which  failure  occurs  outside  the  area  A.     Then  what  we 
observe  is  one  of  the  two  events  [Al-1 ,  Al-2] 

[X=x]n[X<Q]   or  [Q=q]0[X>Q]. 

Then  the  likelihood,  L,  of  the  event 

k  n 


F(x)  =  l-exp[-( 


3 


n  [Xi  <  Qi]  [Xi  =  Xi]  n  [Xj>Qj]  [Qj=qj] 
1=1  j=k+l 


(Al.l) 


where  k  <  n  is  the  random  number  of  failures  observed,  has  the  form 


k  n 
L  =  C  n  f(xi)  n  [l-F(qj)] 
i=l  j=k+l 


(A1.2) 


where  f  is  the  density  function,  F  is  the  distribution  function,  and 


k 


n! 


(n-k) !k! 


The  logarithm  of  the  likelihood  function  for  equation  A1.2  is 


k  n 
in  L  =  in  [C  n  f(xi)  E  (l-F(qj))] 
i=l  j=k+l 


and  can  be  rewritten  as 


k 

Jin  L  =  An  C  +  kilna  -  kciAng  +  (a-1)  I  iln(xi-X) 

i=l 

k      x.-X  ^        n        q.-X  ^ 
-  2    (^— )    -    S  (-^) 
i*l        ^  j=k+l  e 

On  differentiating  equation  A1.3  with  respect  to  a,  3,  and  X  in  turn  and  equating 
the  resulting  expressions  to  zero,  one  obtains  the  maximura  likelihood  estimation 
equations  for  the  three  Weibull  parameters  [Al-3].     These  three  equations  are 
nonlinear  in  the  three  parameters  and  can  be  solved  iteratively  by  a  Newton- 
Raphson  method  as  discussed  in  [Al-1]   and  [Al-3]. 
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APPENDIX  II.     DEPENDENCE  OF  ESTIMATED  60- SEC  LOAD  CAPACITY,  P50 (0-008),  UPON  NUMBER 
OF  RING-ON-RING  SPEQMENS  BEING  TESTED 


Estimates  of  the  order  of  magnitude  of  the  error  in  the  estimation  of  the  load 
capacity  p5o(0.008)  were  carried  out  as  follows.     It  was  assumed  that  the 
probability  distribution  of  the  strength  can  be  modeled  by  a  three-parameter 
Weibull  distribution  with  scale,  location,  and  shape  parameters  (corresponding 
to  an  area  A  =  1.56  in'^)  S^  =  118.1  MPa,  \i^=U7 ,7  MPa,  and  m  =  1.27  (see  table 
3.5,  case  3).     It  is  further  assumed  that  the  probability  distribution  of  the 
load  capacity  p^g  is  also  Weibull  with  parameters  to  be  determined. 

The  first  step  in  obtaining  the  desired  estimates  consists  of  generating  from 
this  distribution  by  Monte  Carlo  simulation  a  sample  of  N  values  of  the  glass 
strength  S.  These  N  values  can  be  viewed  as  strengths  that  would  be  measured 
if  N  ring-on-ring  specimens  were  tested. 

The  second  step  is  to  fit  a  three-parameter  Weibull  distribution  to  these  N  values. 

The  third  step  is  to  estimate  from  this  distribution  the  corresponding  60-sec  load 
capacity  p5o(0.008)  as  shown  in  section  2.3. 

These  three  steps  were  carried  out  20  times  for  each  of  the  sample  sizes  N  = 
15,  100,  250,  500,  and  1,000.     The  results  of  the  calculations  are  summarized 
in  table  A2.1,  which  shows  for  each  value  N  the  means  of  20  estimates  p5o( 0.008), 
their  standard  deviation  s[p5o(0.008) ]  ,  and  their  minimum  and  maximum  values, 
min  [p5o(0.008)]  and  max  [(p5o(0.008)] .     Table  A2.1  suggests  the  order  of 
magnitude  of  the  sample  size  needed  to  obtains  estimates  of  P50 (0.008)  with 
various  precisions.     Note  that  these  estimates  are  tentative,  since  (1)  they  do 
not  account  for  the  effect  of  censoring  discussed  in  appendix  Al ,  and  (2)  they 
assume  that  the  Weibull  distribution  parameters  of  the  glass  strength  obtained 
from  the  testing  of  only  89  specimens  are  the  "true"  parameters. 


Table  A2.1     Estimated  Means,  Standard  Deviations,  Minimum  Values,  and  Maximum 

Values  Obtained  from  Sets  of  Estimated  Values  p5o(0.008).  Corresponding 
to  Various  Ring-on-Ring  Test  Sample  Sizes  n. 


15 


100 


250 


500 


1000 


P50(0.008) 


s[p6o(0.008) 


min  [p5o(0.008)] 


max  [p5o(0.008)] 


33.10 


lA.lO 


2.53 


58.24 


34.36 


8.20 


17.54 


48.26 


35.50 


5.16 


26.35 


47.75 


34.98 


4.23 


27.84 


43.45 


34.69 


2.01 


29.53 


37.54 


APPENDIX  III.     COMPUTER  PROGRAM  FOR  ESTIMATING  PROBABILITY  DISTRIBUTION 
OF  LOAD  CAPACITY  p^Q 


C  Computer  program  PSIXTY  is  a  modified  and  corrected  version 
C  of  the  program  listed  in  "Wind  Loading  and  Strength  of 
C  Cladding  Glass"  by  D.A.  Reed  and  E.  Simiu  (NBS  BSS  154). 
C  It  estimates  the  probability  of  failure  of  square  glass 
C  panes  subjected  to  a  constant  sixty  second  pressure.  The 
C  initial  strength  values  for  each  element  on  each  pane  are 
C  simulated  from  Weibull  distribution  input  parameters.  The 
C  stress  at  the  center  of  each  element  is  obtained  as  a  non- 
C  linear  function  of  pressure  by  using  the  Texas  Tech 
C  University  program  referred  to  in  Section  2.3  of  this 
C  report.    The  stress-pressure  relationship  is  summarized 
C  in  subroutine  LF  of  this  program. 
C 

C  To  obtain  failure  statistics,  numerical  experiments  are 

C  conducted.  A  user-defined  number  of  panes  (=N  in  main  program) 

C  are  loaded,  and  a  distribution  of  the  breaking  pressure  is 

C  determined.    The  breaking  pressure  is  different  for  each  pane, 

C  since  each  pane  has  a  random  Weibull  strength  distribution. 

C  The  pane  is  divided  into  square  elemental 

C  areas,  wherein  each  element  has  a  random  Weibull  distribution 
C  strength  assigned  for  each  of  the  20  directions.  Symmetry  is 
C  taken  advantage  of  in  this  program. 
C 

C  This  program  was  developed  for  panes  with  size  48"xM8"x1/8". 
C  However,  with  small  modifications,  it  can  be  used  for  panes 
C  of  any  specified  size. 
C 

C  The  basic  procedure  is: 
C 

C  FOR  EACH  OF    N    PANES  ... 

C  GENERATE  THE  RANDOM  STRENGTHS  FOR  EACH  DIRECTION  IN  EACH 

C  OF  THE  DISCRETE  ELEMENTS  OF  THE  PANE. 

C 

C  FOR  EACH  LOCATION,  FIND  THE  60-SECOND  PRESSURE  THAT  WILL 

C  CAUSE  THE  STRESS=STRENGTH. 

C 

C  FIND  THE  MINIMUM  60-SECOND  PRESSURE  FOR  THE  PANE. 

C 

C  FIT  THE  N  MINIMUM  60-SECOND  PRESSURES  TO  A  WEIBULL  DISTRIBUTION. 
C 

C  INPUT:    N=Number  of  panes  to  use. 

C  SSHAPE,SSCALE,SLOC=  Input  Weibull  strength  distribution 

C  parameters. 

C  SEEDsRandom  seed  for  random  #  generator. 

C 

C  OUTPUT:  60-second  breaking  pressure  (psf )  parameters  of  the 
C  fitted  Weibull  distribution. 

C  60-second  breaking  pressure  statistics  (mean,S. D. ,C.  0. V. ) 

C 

C  Subroutines  called: 
C 

C     Internal:  F 


C  WEIB 

C  PSIXTY 

C  FDMIN 

C  DIR 


C  LF 
C 

C  External:  GGUBS  -  an  IMSL  routine  to  generate  uniform  random  #'3 
C  SRTAD  -  an  IMSL  routine  to  sort  inplaoe  the  elements 

C  of  a  vector. 

C 

REAL  MEAN 


DOUBLE  PRECISION  SEED, SEEDS 
COMMON  /P60PAR/  PSHAPE, PSCALE, PLOC 
READ  (5,100)  N 

READ  (5,110)  SSHAPE,SSCALE,SLOC,AREA 
READ  (5,120)  SEED 


100  FORMAT  (110) 
110  FORMAT  (4F10.0) 
120  FORMAT  (D20.0) 


SEEDS  =  SEED 


WRITE  (6,130)  SSH APE, SSCALE,SLOCy AREA, N, SEEDS 
130  FORMAT  (//, 

146H  COMPUTER      PROGRAM      PSIXT  Y,//, 

2  37H  GLASS  STRENGTH  WEIBULL  DISTRIBUTION:,/. 

3  21H  SHAPE  PARAMETER        =,  F8.3,/, 

3  21H  SCALE  PARAMETER        =,  F7.2,/, 

4  21H  LOCATION  PARAMETER  =,  F7.2,//, 

5  21H  TEST  SPECIMEN  AREA  =,  F8.3,//, 

6  21H  NUMBER  OF  PANELS      =,  15,//, 

7  21H  RANDOM  NUMBER  SEED  =,  D20. 10) 


CALL  PSIXTY 

+  ( N , AREA , SSHAPE , SSCALE , SLOC , PSHAPE , PSCALE , PLOC , MEAN , STDDE V , COV , 
+  SEED) 


F8  =  F(0.008) 
FM  =  F(0.5) 


WRITE  (6,140)  PSHAPE, PSCALE, PLOC, ME AN, STDDEV, COV, F8.FM 
140  FORMAT  (//, 

1  19H  PSIXTY  PARAMETERS:,/, 

2  21H  SHAPE  PARAMETER        =,  F10.3,/, 

3  21H  SCALE  PARAMETER        =,  F  9.2,/. 

4  21H  LOCATION  PARAMETER  =,  F  9.2,//, 

5  21H  MEAN  =,  F  9.2,/, 

6  21H  STANDARD  DEVIATION  =,  F10.3,/, 

7  21H  COV  =,  FIO.3,//, 

8  21H  F(0.008)  =,  F  9.2,/. 

9  21H  F(0.5)     (MEDIAN)      =,  F  9.2) 


END 

FUNCTION  F(Z) 


COMMON  /P60PAR/  PSHAPE, PSCALE, PLOC 


F  =  PLOC  +  PSCALE«(-ALOG(1.0  -  Z) )••( 1. 0/PSHAPE) 


RETURN 
END 

SUBROUTINE  WEIB(N,M,MR,SSHP,SSCL,SL0C,E1 ,E2,E3,T,S1 ,S2,S3) 

C  PROGRAM  WEIBULL— COMPUTES  THREE  PARAMETER  WEIBULL 

C  MAXIMUM  LIKELIHOOD  ESTIMATES  FOR  BOTH  LEFT  AND 

C  RIGHT  SINGLE  CENSORING 

C  DATE:    JUNE  31,1981 

C  DATE  PROGRAM  CHECKED  OUT:    JULY  9,1981 

C  PROGRAMMER:    JONATHAN  W.  MARTIN 

C  REFERENCE:    HARTER,H.  L. ;MOORE, A. H.     MAXIMUM  LIKELIHOOD 

C  ESTIMATION  OF  THE  PARAMETERS  OF  GAMMA  AND  WEIBULL 

C  POPULATIONS  FROM  COMPLETE  AND  CENSORED  SAMPLES. 

C  TECHNOMETRICS  7: 639-643; 1965.    ERRATA, 9: 195; 1967. 

C  NOTATION  FOR  INPUT  DATA: 

C  N    =    SAMPLE  SIZE  (BEFORE  CENSORING),  N=1000  OR  LESS 

C  AS  DIMENSIONED 

C  SSCL=    0  IF  THE  SCALE  PARAMETER  THETA  IS  KNOWN 

C  SSCL  =  1  IF  THE  SCALE  PARAMETER  THETA  IS  UNKNOWN 

C  SSHP  =  0  IF  THE  SHAPE  PARAMETER  K  IS  KNOWN 

C  SSHP  =  1  IF  THE  SHAPE  PARAMETER  K  IS  UNKNOWN 

C  SLOC  =  0  IF  THE  LOCATION  PARAMETER  C  IS  KNOWN 

C  SLOC  =  1  IF  THE  LOCATION  PARAMETER  C  IS  UNKNOWN 

C  T(I)  =  I-TH  ORDER  STATISTIC  OF  SAMPLE  (I=1,N) 

C  (SUBSTITUTE  BLANK  CARDS  FOR  UNKNOWN  CENSORED 

C  OBSERVATIONS) 

C  M      =    NUMBER  OF  OBSERVATIONS  REMAINING  AFTER 

C  CENSORING  N-M  FROM  ABOVE 

C  C(1)=INITIAL  ESTIMATE  (OR  KNOWN  VALUE)  OF 

C  THETA(I)  =  INITIAL  ESTIMATE (OR  KNOWN  VALUE)  OF  THETA 

C  EK(1)  =  INITIAL  ESTIMATE  (OR  KNOWN  VALUE  )  OF  K 

C  MR  =  NUMBER  OF  OBSERVATIONS  CENSORED  FROM  BELOW, 

C  NORMALLY  ZERO  INITIALLY 

C  NOTATION  FOR  OUTPUT  DATA: 

C  N, SSCL, SSHP, SLOC, M,C(1),THETA(1),EK(1)— SAME  AS  FOR  INPUT 

C  C(J)  =  ESIMATE  FOR  LOCATION  PARAMETER  AFTER  J-1 

C  ITERATIONS  (OR  KNOWN  VALUE  )  OF 

C  THETA(J)  =  ESTIMATE  OF  SCALE  PARAMETER  AFTER  J-1 

C  ITERATIONS  (OR  KNOWN  VALUE)  OF  THETA 

C  EK(J)  =  ESTIMATE  OF  SHAPE  PARAMETER  AFTER  J-1 

C  ITERATIONS  (OR  KNOWN  VALUE)  OF  K 

C  (MAXIMUM  VALUE  OF  J  AS  PRESENTLY  DIMENSIONED 

C  IS  550) 

C  EL  =  NATURAL  LOGARITHM  OF  LIKELIHOOD  FOR  C(J) ,THETA(J) , 

C  EK(J) 


DIMENSION  THETA(5550),EK(5550) 
DOUBLE  PRECISION  C ( 5550) ,T( 1000) ,SLK 

DIMENSION  X(600),Y(600) 

REAL  E1,E2,E3 


C  INTEGER  PSCL(2),PSHP(2),PL0C(2) 

C  DATA  PSCL/'  YES','  NO  '/, 

C  1  PSHP/'  YES','  NO  '/, 

C  2         PLOC/'  YES','  NO  '/ 

C  WRITE(6,«)N,M,MR,SSHP,SSCL,SL0C 

C  WRITE(6,«)E1,E2,E3 

C    4  READ(5,1)N,M,MR,SSHP,SSCL,SL0C 

C       IF(N)  66,66,77 

C    77  READ(5,5)  EK( 1 ) ,THETA( 1 ) ,C( 1 ) 

EK(1)=E1 

THETA(1)=E2 

C(1)=E3 
C      5  FORMAT (3F5.0) 
C  READ(5,121)(DATAS(J),J=1,18) 

ISCL  =  SSCL  +  1 

ISHP  =  SSHP  +  1 

ILOC  =  SLOC  +  1 

EN  =  N 
C         DO  2  1=1, N 

IF  (M)  64,64,32 
32  EM  =M 
31  ELNM  =  0. 

JCOUNT  =  0 

EMR  =  MR 

MRP  =  MR  +  1 

NM  =  N  -  M  +  1 

DO  34  I=NM,N 

EI  =  I 

34  ELNM  =  ELNM  +  ALOG(EI) 
IF(MR)  66,35,74 

74  DO  75  1=1, MR 
EI  =  I 

75  ELNM  =  ELNM  -  ALOG(EI) 

35  DO  30  J=1,550 
JCOUNT  =  JCOUNT  +  1 
IF(J-I)  66.25,37 

37  JJ  =  J-1 
SK  =  0. 
SL  =  0. 

DO  6  I  =  MRP,M 

6  SK  =  SK  +  (T(I)-C(JJ))«*EK(JJ) 
IF(SSCL)  7,7,8 

7  THETA(J)  =  THETA(JJ) 
GO  TO  9 

8  IF(MR)  66,19,20 

19  THETA(J)  =  ((SK  +(EN-EM)»(T(M)-C(JJ))»«EK(JJ))/EM) 
1»»(1./EK(JJ)) 

GO  TO  9 

20  X(1)  =  THETA(JJ) 
LS  =  0 

DO  21  L=1,55 
LL  =  L  -  1 
LP  =  L  +  1 
X(LP)  =  X(L) 

ZRK  =  ((T(MRP)-C(JJ))/X(L))««EK(JJ) 


Y(L)  =  -EK(JJ)»(EM-EMR)/X(L)  +  EK(JJ)»SK/X(L)»«(EK( JJ)  + 
1+EK(JJ)«(EN-EM)«(T(M)-C(JJ))««EK(JJ)/X(L)»»(EK(JJ)+1. ) 
2-EMR«EK ( J J ) »ZRK»EXP ( -ZRK ) / ( X ( L) • ( 1 . -EXP ( -ZRK ) ) ) 

IF(Y(L))  53,73,54 

53  LS  =  LS  -1 

IF(LS  +  L)  58,55,58 

54  LS  =  LS  +  1 

IF(LS  -  L)  58,56,58 

55  X(LP)  =  .5»X(L) 
GO  TO  61 

56  X(LP)  =  1.5»X(L) 
GO  TO  61 

58  IF(Y(L)»Y(LL))  60,73,59 

59  LL  =  LL-1 
GO  TO  58 

60  X(LP)  =  X(L)  +  Y(L)«(X(L)-X(LL))/(Y(LL)-Y(L)) 

61  IF(ABS(X(LP)-X(L))-1.E-4)  73.73,21 
21  CONTINUE 

73  THETA(J)  =  X(LP) 
9  EK(J)  =  EK(JJ) 
IF(SSHP)  12,12,11 

11  DO  17  I  =  MRP,M 

17  SL  =  SL  +  DLOG(T(I)-C(JJ)) 
X(1)=EK(J) 

LS=0 

DO  51  L=1,55 

SLK  =  0. 

DO  18  I=MRP,M 

18  SLKsSLK  +(DLOG(T(I)-C(JJ))-ALOG(THETA(J)))»(T(I)-C(JJ)) 
1»«X(L) 

LL  =  L-1 
LP  =  L+1 
X(LP)=X(L) 

ZRK=((T(MRP)-C(JJ))/THETA(J))»«X(L) 

Y ( L ) = ( EM-EMR ) « ( 1 . /X ( L ) -ALOG ( THETA ( J ) ) ) +SL-SLK/THETA ( J ) 
1«»X(L)+(EN-EM)»(AL0G(THETA(J))-DL0G(T(M)-C(JJ))) 
2»(T(M)-C(JJ))»»X(L)/THETA(J)««X(L) 
3+EMR«ZRK» ( ALOG ( ZRK ) /X ( L ) ) •EXP ( -ZRK ) / ( 1 . -EXP ( -ZRK ) ) 

IF(Y(L))  43,52,44 

43  LS  =  LS  -1 
IF(LS+L)  47,45,47 

44  LS  =  LS  +  1 
IF(LS-L)  47,46,47 

45  X(LP)  =  .5»X(L) 
GO  TO  50 

46  X(LP)  =  1.5«X(L) 
GO  TO  50 

47  IF(Y(L)«Y(LL))  49,52,48 

48  LL  =  LL  -  1 
GO  TO  47 

49  X(LP)  =  X(L)+Y(L)«(X(L)-X(LL))/(Y(LL)-Y(L)) 

50  IF(ABS(X(LP)-X(L))-1.E-4)  52.52,51 

51  CONTINUE 

52  EK(J)  =  X(LP) 

12  C(J)  =  C(JJ) 
IF(SLOC)  25,25,14 


14  IF(1.-EK(J))  16,78,78 
78  IF(SSCL+SSHP)  57,57,16 
16  X(1)  =  C(J) 

LS  =  0 

DO  23  L=1,55 

SKI  =  0. 

SR  =  0. 

DO  15  I=MRP,M 

SKI  =  SKI  +(T(I)-X(L))»«(EK(J)-1.  ) 

15  SR  =  SR+  1./(T(I)-X(L)) 
LL  =  L-1 

LP  =  L  +  1 
X(LP)  =  X(L) 

ZRK  =  ((T(MRP)  -X(L))/THETA(J))»«EK(J) 

Y(L)  =  (1.-EK(J))«SR+EK(J)»(SKU(EN-EM)«(T(M)-X(L)) 

1»«(EK(J)-1. ))/THETA(J)««EK(J)-EMR»EK(J)»ZRK» 

2EXP(-ZRK)/((T(MRP)-X(L))«(1.-EXP(-ZRK))) 
IF(Y(L))  39,24,40 

39  LS  =  LS-1 
IF(LS+L)  70,41,70 

40  LS  =  LS+1 
IF(LS-L)  70,42,70 

41  X(LP)  =  .5»X(L) 
GO  TO  22 

42  X(LP)  =  .5»X(L)+.5»T(1) 
GO  TO  22 

70  IF(Y(L)«Y(LL))  72,24,71 

71  LL  =  LL  -  1 
GO  TO  70 

72  X(LP)  =  X(L)+Y(L)»(X(L)-X(LL))/(Y(LL)-Y(L)) 

22  IF(ABS(X(LP)-X(L))-l.E-4)  24.24,23 

23  CONTINUE 

24  C(J)  =  X(LP) 
GO  TO  25 

57  C(J)  =  T(1) 

25  IF(MR)  66.38,69 
38  DO  63  I  =  1,M 

IF(C(J)+1.E-4-T(I))  68,67,67 

67  MR  =  MR  +  1 
WRITE(6,201) 

201  FORMAT (»   ',»  LINE  200  IN  PROGRAM— INITIAL  LOCATION  PARAMETER', 
1 'ESTIMATE  GREATER'/ 'THAN  FIRST  OBSERVED  FAILURE  TIME') 

63  C(1)  =  T(1) 

68  IF(MR)  66,69,31 

69  SK  =  0. 
SL  =  0. 

DO  36  I=MRP,M 
SK  =  SK+(T(I)-C(J))»«EK(J) 
36  SL  =  SL+DLOG(T(I)-C(J)) 

ZRK=((T(MRP)-C(J))/THETA(J))»«EK(J) 

EL=ELNM+ ( EM-EMR ) • ( ALOG ( EK ( J ) ) -EK ( J ) » ALOG ( THETA ( J ) ) ) + 
1(EK(J)-1.)»SL-(SK+(EN-EM)«(T(M)-C(J))«»EK(J))/ 
2 ( THETA ( J ) »»EK ( J ) ) +EMR* ALOG ( 1 . -EXP ( -ZRK ) ) 

IF(J-3)  30.27,27 
27  IF(ABS(C(J)-C(JJ))-1.E-4)  28.28,30 


28  IF(ABS(THETA(J)-THETA(JJ))-1.E-4)  29.29,30 

29  IF(ABS(EK(J)-EK(JJ))-1.E-4)  99.99,30 

30  CONTINUE 


C  OUTPUT 

99  CONTINUE 

S1=THETA(JC0UNT) 
S2=EK(JC0UNT) 

S3=C(JC0UNT) 
GO  TO  66 
64  CONTINUE  . 
66  RETURN 
END 

SUBROUTINE  PSIXTY ( NPANEL , AREA , C 0 , SO , SU , MP , PEST , SV 
1  .XMEAN,ZSAM,COV,DSEED) 

C 

C    (SEE  DATA  STATEMENT)  v 
C  SIDEA  =  LENGTH  OF  THE  SQUARE 

C  THE  LENGTH  OF  THE  SIDE  OF  THE  SQUARE  IS  INPUT  IN  INCHES 

C  TH  =  THICKNESS  OF  THE  PLATE 

C  THICKNESS  INPUT  IN  INCHES 

C  MODULUS  OF  ELASTICITY  (E)  INPUT  IN  PSI 

C 

C  INPUT: 

C    NPANELS  =  #  PANELS  TO  'BREAK*  TO  FIT  WEIBULL  DISTR.  TO 

C  SO,CO,SU  =  WEIBULL  STRENGTH  PARAMETERS  FOR  AN  ELEMENTAL  AREA  = 

C 

C  OUTPUT: 

C  XMEAN  =  MEAN  P60 
C  ZSAM    =  S.D.  OF  P60 
C  COV      =  C.O.V.  OF  P60 

C  PEST    =  SCALE  PARAMETER  OF  WEIBULL  MLE  FIT  OF  P60  FOR  NPANELS 
C  MP        =  GAMMA    "  n  n 

C  SV        =  LOCATION  PARAMETER    "  " 
C 

C  PROCEDURE: 

C  FOR  EACH  OF  NPANEL  PANES... 

C  GENERATE  216  R.V.  STRENGTHS  (FROM  UNDERLYING  WEIBULL 

C  DISTR. )  36  LOCATIONS  W/  EACH  6  DIRECTI0NS=216 

C 

C  FOR  EACH  LOCATION,  IF  SIGMA=STRENGTH,  FIND  CORRESPONDING 

C  LOAD(P  60)  WHICH  GIVES  THAT  STRESS ( STRENGTH ) ,  I.E.  THAT 

C  LOAD  WHICH  WILL  BREAK  IT  FOR  THAT  LOCATIONS  STRENGTH. 

C 

C  FIND  THE  MINIMUM  LOAD  (P60)  FOR  THIS  PANEL. 

C 

C    FIT  THE  MINIMUM  LOADS  (P  60)  FOR  THE  NPANEL  PANES  TO  A  WEIBULL 

C  DISTRIBUTION. 

C 

C  NOTE  THAT  SYMMETRY  IS  USED  IN  THIS  PROCEDURE,  I.E.  THE  PANE  IS 
C  DIVIDED  INTO  FORU  QUADRANTS  (BY  SYMMETRY)  AND  THIS  QUADRANT  IS 
C    DIVIDED  BY  ITS  AXIS  OF  SYMMETRY  (THE  DIAGONAL)  IN  TWO. 
C 

C  THE  NUMBERING  IN  SUBROUTINE  DIR  AND  LF  USES  THE  UPPER  HALF  OF 
C  LOWER  LEFT  QUADRANT.     THE  CENTROID  OF  THE  PANE  IS  NUMBERED 


C  NODE  #1.     THE  24"  SQUARE  QUADRANT  OF  THE  MS"  SQUARE  PANE  IS 

C  DIVIDED  INTO  8X8  3  INCH  SQUARES.  NODE  #  1  IS  THE  CENTROID  OF 

C  THE  ENTIRE  PANE  AND  IS  THE  CENTROID  OF  ELEMENT  #1.  NODE  #  2  IS 

C  3  INCHES  LEFT  OF  NODE  #1  AND  IS  THE  CENTROID  OF  ELEMENT  #2. 

C  ETC.  TILL  NODE  #  8  IS  21  INCHES  LEFT  OF  NODE  #1  AND  IS  THE 

C  CENTROID  OF  ELEMENT  #8.  AGAIN,  USING  SYMMETRY,  NODE  #9  IS 

C  DIRECTLY  UNDER  N0DE#2  AND  IS  THE  CENTROID  OF  ELEMENT  #9. 

C  NODE  #10  IS  3  INCHES  LEFT  OF  NODE  #9  ETC  UNTIL  NODE  #  15 

C  WHICH  IS  18  INCHES  LEFT  OF  NODE  #9.  THEN  NODE  #  16  IS  DIRECTLY 

C  UNDER  NODE  #10.  ETC.  THEN  NODE  #22  IS  DIRECTLY  UNDER  NODE* 17. 

C  ETC.  UNTIL  THE  LOWEST  ROW  WHICH  HAS  ONLY  NODE  #36. 

C  NOTE  THAT  EACH  ELEMENT  IS  A  THREE  INCH  SQUARE  (THIS  IS 

C  NECESSARY  SINCE  THE  R.V.  STRENGTH  ASSUMES  AN  EQUAL  ELEMENTAL 

C  AREA.  THUS,  ELEMENT  #1  HAS  NO  COUNTERPART  IN  THE  OTHER  THREE 

C  QUADRANTS.  LIKEWISE  ELEMENT'S  #2-8  EACH  REPRESENT  FOUR  OTHER 

C  ELEMENTS  IN  THE  PANE  (ALSO  WITH  THE  ELEMENTS  ON  THE  DIAGONAL) 

C  THE  REMAINING  ELEMENTS  ARE  REPLICATES  OF  8  OTHER  ELEMENTS  IN 

C  THE  PANE. 

C 

DOUBLE  PRECISION  SIG(36 ,6) ,SI (36 ,6) ,PEQ( 36 ,6) , A1 (216) , 
•RSIG ( 36 , 6 ) , S2 , BP , M , RN , RY , RNM , KIC , ROOT , COEFS , COEFP , 
•SCALE , SMIN , SMAX , RF , DSEED , E . MIN ( 1 000 ) 

REAL  F(216),ZSAM,C0V,XMEAN,MP 

INTEGER  IPOINT,NUPPER 

DIMENSION  ICHEK(216,2) 

EQUIVALENCE  (SI(36,6) ,PEQ(36,6) ) 

COMMON  //IPLACE(216,2),IC0UNT 

C 

DATA  KIC , A , RN , RY/ 0. 75D0 , 1 . 08 , 1 9. 6 9D0 , 1 . 25D0/ 

DATA  TH,SIDEA,PR/0. 125,48.0,0.21/ 

DATA  E/1D07/ 

DATA  RF,IPOINT/1.0D0,1/ 

C 

C  CALCULATIONS 
C 

M=1/C0 
RNM=RN-2. 
ROOTS 1./RN 

C 

C  1/B  IS  A  PARAMETER  USED  IN  EQN.   19  OF  THE  BSS  REPORT 
C  RY=1.25;  A=1.08;  N=19.69  FROM  [5]  OF  BSS  154 
C  SEE  ABOVE  DATA  STATEMENT 
C 

BP=(RNM«A«RY»RY»(KIC«»RNM) )/2. 

C 

C  FLEX=FLEXURAL  RIGIDITY 
C 

FLEX=(  E»TH»TH«TH  )/(  12. »( 1. -PR«PR  )) 

C 

S2=SIDEA«SIDEA 

C 

C  CF.  EQN  (9)  OF  BSS  154 

C  THE  AREA  FOR  WHICH  WE  ARE  SIMULATING  A  STRENGTH  IS  THE 
C  ELEMENT  SIZE,  WHICH  AS  EXPLAINED  ABOVE  IS  1/16  OF  SIDE A 
C  THUS,  AREA=  (SIDEA/16)»»2 

C  HOWEVER,  TAKING  DIRECTION  INTO  ACCOUNT  AS  EXPLAINED  BELOW, 


C  THIS  AREA  IS  DIVIDED  BY  10.0 
C 

SELEM=(SIDEA/16.  )«»2 

SNEW=SO»  (AREA/(SELEM/10. ))*«M 

C0EFS=(S2«TH)/FLEX 

COEFP= ( S2»S2 ) / (FLEX«TH ) 

C 

C  THESE  COEFFICIENTS  ARE  USED  FOR  NONDIMNSIONALIZING 

C  PRESSURES  AND  STRESSES. 

C 

SCALES ( 1 . / ( 60. »BP ) ) «*ROOT 
R00T=RNM«R00T 

C 

C    OBTAIN  THE  DIRECTIONAL  MULTIPLIERS  FROM  [9]  OF  BSS  154 
C 

CALL  DIR(RSIG) 

C 

C  THEN  CALCULATE  THE  STRESSES  IN  TERMS  OF  THE  STRENGTHS... 
C 

C  IF  THE  STRESSES  ARE  LESS  THAN  OR  EQUAL  TO 
C  ZERO,  THEN  OMIT  THESE  FROM  CONSIDERATION... 
C 

C  LOGICAL  UNIT  9  CONTAINS  A  FILE  WITH  PSEUDO-RANDOM 

C  NUMBERS  WHICH  FOLLOW  A  UNIFORM  DISTRIBUTION  ON 

C  (0,1).  AN  ALTERNATIVE  METHOD  WOULD  BE  TO  GENERATE  THE  NUMBERS 

C  DIRECTLY  FROM  A  SUBROUTINE. 

C 

C  THEN  GENERATE  INITIAL  STRENGTH  VALUES 
C 

DO  926KK=1,NPANEL 
1=0 

DO  12J=1,36 
IF(J.EQ.  1)  THEN 
Nsl 

ELSEIF(J.EQ.9  .OR.  J.EQ.  16  .OR.  J.EQ.22  .OR.  J.EQ.27  .OR. 
+    J.EQ.  31  .OR.  J.EQ.  34  .OR.  J.  EQ.  36  .OR.   ( J. GE.  2. AND.  J. LE. 8) )THEN 

N=4 
ELSE 

N=8 
ENDIF 

C 

C  SINCE  THE  STRENGTH  IS  RANDOMLY  SIMULATED  W/IN  EACH  ELEMENT 
C  WE  CAN  ARBITRARILY  ORIENT  THE  COORD  SYSTEM,  S.T.  THE  MAXIMUM 
C  PRINCIPAL  STRESS  OCCURS  §  ALPHA=0 

C  DIVIDING  THE  CIRCLE  INTO  20  =  18  DEGREE  ARCS,  THE  FIRST  ARC 

C  IS  REPEATED  TWICE,  @  ARC#1  AND  180  DEGREES  AROUND  THE  CIRCLE. 

C  SIMILARLY,  ARC#6,  CORRESPONDING  TO  THE  MINOR  PRINCIPAL  STRESS 

C  ALSO  OCCURS  TWICE.  THE  STRESSES  IN  THE  ARCS  BETWEEN  ARC  #  1  AND 

C  ARC  #  6  OCCUR  FOUR  TIMES  (BY  SYMMETRY).  THUS  STRESS  @  ARC#1 

C  OCCURS  1/10;  STRESS  e  ARC#6  OVVURS  1/10;  REMAINING  ARCS 

C  STRESS  OCCURS  1/5  EACH;   1/10+1/10+4( 1/5)=1. 0 

C  USING  THE  ANALOGY  OF  EQN  (11)  IN  BSS  154 

C  THE  AREA  CORRESPONDING  TO  SO  IS  3  SQUARE  INCHES  DIVIDED 

C  BY  10.  ALSO,  THE  STRESSES  CORRESPONDING  TO  ARC#2-5, 

C  ARE  SIMULATED  TWICE,  BECAUSE  THEY  OCCUR  TWICE  AS  OFTEN  AS 

C  THOSE  CORRESPONDING  TO  THE  MAJOR  AND  MINOR  PRINCIPAL  STRESSES 


C  THUS.... 

C 

DO  12  K=1,6 

IF  (K.NE.  1  .AND.  K.NE.6)  N=N«2 

C 

C  GGUBS  is  an  IMSL  routine  to  generate  N  random  uniform  #'s 
C 

CALL  GGUBS (DSEED,N,F) 
CALL  FDMIN(F,N,FMIN) 

IF  (  FMIN  .EQ.   0  )  THEN 

FMIN=FMIN+0. 00001 

ENDIF 

C  SI=  STRENGTH  (UNITS  OF  MPA) 

C  SIG=STRESS  FACTOR  (SF)  (UNITS  OF  PSI)  =  F(SIGMA  60) 
C  SIG  FROM  EQN.   (13)  &  (19)      0.00689  PSI  PER  MPA 
C 

C  FOLLOWING  LINE  CHANGED  FROM  FMIN      TO  (1.-FMIN) 

C  ORIG.  WAS  OK  SINCE  ONLY  ONE  RANDOM  #0-1     &    A  RANDOM  UNIFORM  # 

C  BETWEEN  0-1  IS  EQUIVALENT  TO  RANDOM  1-  RANDOM  # 

C  C 

SI(J,K)=SNEW»((-AL0G((1.-FMIN)))»»M)  +SU 
SIG(J,K)=SCALE»(  SI(J,K)»»R00T  ) 
SIG(J,K)=  C0EFS«(  SIG(J,K)/(  . 00689»RSIG( J.K) ) ) 
IF(  SIG(J,K)  .GT.  0)THEN 
1=1+1 
C  J==>  POSITION 
C  K==>  DIRECTION 

ICHEK(I,1)=J 
ICHEK(I,2)=K 
ENDIF 
12  CONTINUE 
ICOUNT=I 

C 

C  SORT  THE  REMAINING  STRESSES.... 
C 

DO  160  I=1,IC0UNT 
J=ICHEK(I,1) 
K=ICHEK(I,2) 
A1(I)=SIG(J,K) 
160  CONTINUE 

C 

IC=ICOUNT 

Cx 

C  SRTAD  is  an  IMSL  routine  to  sort  inplace  the  IC  elements  of  vector 

C    A1  min  to  max. 

C 

CALL  SRTAD(A1,1,IC) 
SMIN=A1(1) 
SMAX=A1(IC0UNT) 
SMAX=RF»SMAX 

C 

C  OMIT  FROM  CONSIDERATION  THE  FOLLOWING  STRESS  VALUES... 
C 

GO  TO  (10,10,30)IPOINT 

C 

C  FOLLOWING  LOOP  IS  UNNECESSARY  IF  RF=0.0! 


10  12=0 

DO  161  I=1,IC0UNT 

J=ICHEK(I,1) 

K=ICHEK(I,2) 
IF(  SIG(J,K)  .LE.  SMAX  )THEN 
12=12+1 

IPLACE(I2,1)=J 
IPLACE(I2,2)=K 
ENDIF 

161  CONTINUE 

C 

IC0UNT=I2 

C 

GO  TO  15  - 

C 

30  CONTINUE 

DO  162I=1,IC0UNT 
J=ICHEK(I,1) 
K=ICHEK(I,2) 

IF(  SIG(J.K)  .EQ.  SMIN  )THEN 

IPLACE(1,1)=J 

IPLACE(1,2)=K 
ENDIF 

162  CONTINUE 

C 

C  CALCULATE  LOAD  FACTOR,  LF  ...  PEQ=LF=F(SF). 
C 

15         CALL  LF(SIG,IPOINT,PEQ) 

C 

GO  TO  (11,11,31)IP0INT 

C 

11  JK=0 

DO  17  IV=1,IC0UNT 

J=IPLACE(IV,1) 

K=IPLACE(IV,2) 

JK=JK+1 
IF(  PEQ(J,K)  .LE.  0  )THEN 
PEQ(J,K)=1.0D20 
ENDIF 

C  A1=LF/C0EFP  =  PRESSURE    (UNITS  OF  ?) 
A1(JK)=PEQ(J,K)/C0EFP 
17  CONTINUE 

C 

C  SORT  VALUES  .... 
C 

CALL  SRTAD ( A 1 , 1 , ICOUNT ) 
MIN(KK)=144.»A1(1) 
GO  TO  926 

C 

C  IF  ONLY  THE  MINIMUM  STRESS  IS  CHECKED, 

C  CALCULATE  ONLY  ONE  VALUE  OF  PEQ  

C 

31  J=IPLACE(1,1) 
K=IPLACE(1,2) 

MIN(KK)=144.»PEQ(J,K)/C0EFP 

C 


926  CONTINUE 

C 

C         UNITS  CHANGED  TO  PSF  FOR  CALCULATIONS 

C  MIN=P60(PSF)MIN  OF  PANEL'S  P60 

C 

CALL  SRTAD (MIN , 1 , NPANEL ) 

C 

SUM1=0.0 
SUM2=0. 0 

DO  30031=1, NPANEL 

3003  SUM1=SUM1+MIN(I) 
XMEAN=SUM1 /NPANEL 
DO  30041=1, NPANEL 

3004  SUM2=SUM2+(  MIN (I)  -  XMEAN  )«*2 
ZSAM=SQRT(SUM2/(  NPANEL  -  1  )) 
COV=ZSAM/XMEAN 

C 

N=NPANEL 

C  SSHP,SSCL,SL0C=1  ===>SHAPE, SCALE, LOCATION  PARAMETER  UNKNOWN 

C  ElsINITIAL  ESTIMATE  OF  SHAPE  PARAMETER 

C  E2=  "  n  n  SCALE  " 

C  E3=  "  "  "  LOCATION  " 

C 

SSHP=1. 
SSCL=1. 
SL0C=1. 
E1=2.3 

E2=:10. 

E3=MIN(1)  -  2.0 

IF(E3  .LE,  0.0)  E3=0.001 

IF(  SU  .LE.  0.1)  THEN 

E3  =  0.001    @  SLOC  =  0.0 

ENDIF 
311  MR=0 
C  WRITE (6, 305) 

305    FORMAT ( IX, 3H305) 

C 

CALL  WEIB ( N , N , MR , SSHP , SSCL , SLOC , E 1 , E2 , E3 . MIN . PEST , MP , SV ) 

C 
C 

RETURN 
END 

SUBROUTINE  FDMIN(F,N,FMIN) 
DIMENSION  F(N) 

C 

FMIN=1.E+10 
DO  10  1=1, N 
FMIN=AMIN 1 ( FMIN , F ( I ) ) 
10  CONTINUE 
RETURN 
END 

C 

SUBROUTINE  DIR(RSIG) 

DOUBLE  PRECISION  PH,RSIG(36 ,6) ,C(6) ,SII(6) ,S2S1 (36) 
DATA  S2S1/1.0D0,.98D0,.94D0,.89D0,.79D0,.68D0..52D0,.  12D0,.97D0. 
».93D0,.86D0,.77D0,.65D0,.47D0,.05D0,.89D0,.82D0,.71D0..57D0,.36D0, 


09D0, 0. 73D0, 0. 61BO,iO. 
». 21D0,-. 24D0,, 47D0,. 28D0, , O2D0,-. 38D0, » 83D0,-o 17D0,-«  52D0, 
«-0.39D0,-0.66D0,-.8lD0/ 
DO  5K=1,6 

PH=(K-1)«3. 14159/10. 
C (K ) =DCOS ( PH ) »DCOS ( PH ) 
5  SII(K)=DSIN(PH)«DSIN(PH) 
DO  10J=1,36 
DO  10K=1,6 
RSIG(J,K)=C(K)  +  S2S1(J)»SII(K) 
10  CONTINUE 
RETURN 

END  \  ^ 

C 

SUBROUTINE  LF(SIG, IPOINT, PEQ) 

DOUBLE  PRECISION  PEQ(36 ,6 ) ,SIG(36 ,6 ) ,LIF(25) ,L(36 ,25) 
COMMON  //  IPLACE(2l6,2)sIC0UNT 

C 

DATA  ((L(J,I),I=1,24),Jsl,6)/0,6,8,13.7,20.6,27.5,34.0,40.5, 
«46. 6, 53. 0,58. 3, 63. 0,104. 0,129. 0,149, 167, 182, 195. 207, 218, 
•228,310, 359, 407, 481, 0,6. 7, 13. 5, 20. 4, 27. 0,33.7, 40.4, 46.1, 
•52,57.7,63,103,130,151,170,186,199,211,222,233,319,364,410, 
•484.0,6.5,13,19.7,26.3,32.8,39,45,51,56.7,62,104.134,157, 
•177,194,209,223.235,246,333.384.430,497,0,6,12.3,18.5,24.9. 
•30.9,37.2,43,49,54.2,59,103,136,163,186,206,224,240,255, 
•268.370,432,484,550,0,5.5,11.2,16.8,22.6,28.3,33.8,39,45,50, 
•55,90,135,165,192,215.237,257,275,290,420,509,581,661,0, 
•4.7,9.5,14.3,19.2,24.1,28.9,34,39,43,48,90,125.157,186, 
•212,234,257,278,296,457,580,687,792/ 

C 

DATA  ((L(J,I),Is1,24),J=7,11)/0,3.6,7.2,10.8,14.5,l8.2,21.7, 
•25.5,29.3,33,37,71,102,130,156,180,203,224.244,263,426, 
•563,689,811,0,2.0,4,6.1,8.3,10,12,14,17,19,21,40,59,78, 
•94,109,125,139,153,166,279.377,471,563,0,6.7,13.6,20.4, 
•27.  1,33.7,40.0,45.9,52.2,57.5,62.  1,102.6,129.2,149.7,167.7, 
•183,196.4,209.2,220.6,232.4,319.8,375.9,428.4,504.2,0,6.6, 
•13.3,20,26.6,33,39.2,45.  1,51.2,56.6,61.3,102.8,131.9,154.8, 
•174.9,192.3,207.5,221.3,233.7,245.7,337.  1,397.8,454.0,530.8, 
•0,6.3,12.7,19,25.4,31.3,37.4,42.9,49,54.2,59.1,102.6,134.7, 
•161.5,184.7,205.  1,223.3,239.8,254.6,268.3,373.5,441.5, 
•500. 1,574.2/ 

C 

DATA  ((L(J,I),I=1,24),J=12,15)/0,5.8.11.6,17.3,23.2,28.6, 
•34.4,39.6,44.9,50.4,54.9,98.9,134.2,164.5,191.5,215.4, 
•237.3,257.3,275.8,291.8,424.9,516.9,593.1,672.8,0,5, 
•10,15,20,24.7,29.6,34.4,39,43.6,47.9,89.4,125.1,156.9, 
•186,212.5,236.9,259.8,281. 1,300.3,463. 1,587.7,697,800.7, 
•0.3.9,7.9,10,15.5,19.2,23,26.6,30.  1,33.8.37.4,71.2,101.9, 
•130.5,157.1,182,205.1,227.2,247.9,267.3,434.8,573.2,699.7, 
•820.9,0,2.5,4.9,7.3,9.4,11.8,13.9,16,18,2,20.2,22.2,41.9, 
•60.8,78.6,95.7,111.8,126.9,141.6,155.6,169.  1,285.7,385.8, 
•479.8,571.5/ 

C 

DATA  ((L(J,I),I=1,24),J=16,20)/0,6.6.13.3,19.9,26.4,32.7,38.9, 
•44.4,50.4,55.7,60.3,100.7,128.7,151.0,170.4,187,201.9, 
•216.  1,230.2,242.8,349.3,423.7,490.2,572.4.0,6.5,12,8.19.3, 


•25.6,31.6,37.7,^3.  1,48.7,53.9,58.5,99.7,130.8,157.  1,180.5, 
»201. 1,219.9,236.9,252.8,266.6,381.3,463.6,537.5,621.6,0 
•,6.1,12.1,18.1,23.9,29.5,35.2,40.2,45.6,50.5,55,96.8,131.2, 
•161.7,189.4,214.4,237.3,258.3,278.2,295,435.7,534.2,620, 
•705.3,0,5.5,10.8,16.  1,21.2,26.1,31.2,35.7,40.2,44.8,48.9, 
•88.6,123.6,156.  1,186.3,214.4,240.5,264.7,287.8,308.5, 
•481.8,612. 1,724.5,828.9,0,4.5,8.9,13.2,17.3,21.3,25.  1, 
•28.9,32.7,36.2,39.5,72.1,102.4,131.5,159.3,186.  1,210.8, 
•234. 7 , 257. 5 , 278. 9 , 46  0. 4 , 6  05. 2 , 735. 4 , 856. 9/ 

DATA  ( (L( J, I), 1=1, 24), J=21, 24)70,3.3, 6. 2, 9. 3, 12.  1,14.7,17.2, 
•19.5,21.9,24.  1,26.  1,45.6,63.7,81.5,99,116.2,132.2,148.2, 
•163.5,178.5,307.3,413.9,510.9,603.4,0,5.7,6.5,12.9,19.2, 
•25.3,31.2,37,42.4,47.8,52.7,57.3,97.3,127.3,153.0,175.0, 
•195.213,229.244,391.7,493.4,583.5,676.3,0,5.6,6.3,12.4, 
•18.4,24.2,29.8,35.3,40.4,45.6,50.4,54.7,94.6,126.9,156, 
•184,209,233,255,276,444.6,556.9,655.6,750.8,0,5.3,5.8, 
•11.4,16.9,22.2,27.3,32.3,36.9,41.6,45.9,50.2,88.1,121.6, 
•154.184,214,241,268,294,507.1,648.2,769.1,876.5/ 

DATA  ( (L( J, I), 1=1, 24), J=25, 28)70,4.7,5, 9. 9, 14. 6, 19. 1,23.4, 
•27.6,31.5,35.2,39.1,42.4,74.6,103.9,133,161,190,217, 
•244,269,501.7,662.8,802.8,929.4,0,3.6,3.9,7.7,11.3,14.7, 
•17.8,20.9,23.6,26.4,28.9,31.3,52,70.2,88. 1,106,124,141, 
•159,176,344.9,468.2,576.2,674.6,0,5.7,6.2,12.2,18.2,23.9, 
•29.3,34.7,39.7,44.6,49.3,53.6,92.4,124.1,180,205.3,228.2, 
•250. 1,270.8,288.7,436.2,569,687.7,798.2,0,5.6,5.9,11.8, 
•17.3,22.7,27.9,32.9,37.6,42.3,46.6,50.9,88. 1,120.6,181, 
•210.3,237.9,265.2,291.7,316.4,525.4,680.3,811.5,926.77 

DATA  ( (L( J, I), 1=1, 24), J=29, 32)70, 5. 1,5.4,10.7,15.8,20.5,25. 1, 
•29.6,33.8,38.1,41.8,45.5,78.5,82.3,165,193.9,221.8,249.7, 
•277.7,305.1,548.9,739.9,900.6,1041.0,0,4.2,4.5,8.9,12.9, 
•17,20.9,24.4,27.8,31.2,34.2,37,61.8,76.8,119.9,138.6, 
•156.9,175.5,194.2,213.2,395.9,553.2,688.4,806.5,0,5.4, 
•11.6,17.3,22.6,27.8,32.7,37.2,41.9,46.1,50,87,120,151, 
•172.9,210.2,237.9,265.6,292.3,318.1,538.5,700.9,827,923,0, 
•5.4,11.2,16.5,21.5,26.4,31.2,35.5,39.9,43.9,47.9,83,114. 
•144,172.8,201.7,230,258.286.4,314.7,582.5,809.3,998.4, 
•1157.67 

DATA  ((L(J.I),I=1,24),J=33,36)70,4.7,9.9, 14.7, 19.2,23.6,27.8, 

•31.8,35.7,39.2,40,73,99,122,143.9,165.4,186.3,206.9,227.3, 
•248.2,456.5,658.9,843,1001.9,0,5.4,11.2,16.6,21.7,26.6, 
•31.4,36,40.6,44.8,48.8,85.5,118.148.5,178.3,207.8,236.9, 
•265.9,295.323.7,607.3,865.2,1085.7,1264.5,0,5.  1,10.6, 
•15.8,20.7,25.6,30.  1,34.6,39.2,43.3,47.3,83.5,114.5,143.3, 
•170.6,196.7,222.6,247.7,272.6,296.9,538.2,780.4,1015.6, 
•1229.8,0. ,5. 1,10.8,16. 1,21.4,26.4,31.3,36,40.8,45.3,49.7, 
•89.6,124.4,156.7,187.9,217.2,246.4,274.7,302.4,328.9, 
•587.5,846.4,1113.9,1389.17 

DATA  (LIF(I), 1=1, 25)70, 24. 8, 49.7, 74. 5, 99, 124, 149, 174, 199, 
•223,249,497,745,993,1242,1490,1738,1987,2235,2484,4967, 
•7450,9934, 12417, 1000007 


C  APPROXIMATIONS  FOR  THE  LARGEST  VALUES  OF  SIGMA 

C  ARE  TAKEN  FROM  THE  FOLLOWING  REPORT  : 

C  "PROPOSED  METHOD  FOR  DETERMINING  THE  THICKNESS 

C  OF  GLASS  IN  SOLAR  COLLECTOR  PANELS,"  BY 

C  DONALD  MOORE,  JET  PROPULSION  LABORATORY, 

C  CALIFORNIA  INSTITUTE  OF  TECHNOLOGY,  PASADENA, 

C  CALIFORNIA,  MARCH,  1980. 

C 

DO  6000J=1,29 
6000  L(J,25)=1500 

DO  7000J=30,36 
7000  L(J,25)=5000 
C 

C  CALCULATE  THE  LOAD  FACTOR  FOR  EACH  STRESS 
C 

DO  1000  II=1,IC0UNT 
IF(  IPOINT  .EQ.  3  )G0  TO  300 
J=IPLACE(II,1) 
K=IPLACE(II,2) 
SSrSIG(J,K) 
GO  TO  500 
300  JrIPLACE(1,1) 
K=IPLACE(1,2) 
SS=SIG(J,K) 
500    IF(SS  .GT.  L(J,25))  THEN 
PEQ(J,K)=LIF(25) 
GOTO  1000 
ENDIF 

IF(SS.EQ.O)THEN 

PEQ(J,K)=0 

GOTO  1000 
ENDIF 

IF(SS.LT.  L(J,2))  THEN 

PEQ(J.K)=SS«(LIF(2)-LIF(1))/L(J,2) 
GOTO  1000 
ENDIF 

DO  2001=2,25 
IF(SS  .GT.  L(J.I-I)  .AND.  SS  .LT.  L(J,I))  THEN 

R=(SS-L(J,I-1))«(LIF(I)-LIF(I-1))/(L(J,I)-L(J,I-1)) 
1  +  LIF(I-I) 
PEQ(J,K)=R 
GOTO  1000 
ENDIF 

IF(SS  .EQ.  L(J,I))  THEN 
PEQ(J,K)=LIF(I) 
GOTO  1000 
ENDIF 
200  CONTINUE 

IF(  IPOINT  .EQ.  3  )G0  TO  400 
1000  CONTINUE 
C 

400  RETURN 
END 


APPENDIX  IV,       COMPUTER  PROGRAM  FOR  ESTIMATING  PARAMETERS  OF  GLASS 
STRENGTH  DISTRIBUTION 


*BATCH 


c 

C  PROGRAM  CDMPUTES  THE  MAXIMUM  LIKELIHOOD  ESTIMATES 

C  FOR  THE  TWO  OR  THREE  PARAMTER  WE  I  BULL  AND 

C  FOR  THE  TWO  PARAMETER  LOGNORMAL  FOR  PROGRESSIVE  CENSORED 

C  SAMPLES 

C 

C  OUTPUT   INCLUDES  PARAMETRIC  ESTIMATES,    NUMBER  OF  SAMPLES 

C  AND  NUMBER  OF  FAILED  SAMPLES,      FOR  THE  PURPOSES  OF  THIS 

C  PROJECT,    ONLY   THE   ESTIMATES  OF   THE  WE  I  BULL  DISTRIBUTION 

C  ARE  USED. 
C 

C  SUBROUTINE  HAMLET   IS  USED  TO  OBTAIN  TWO  PARAMETERS 

C  WHILE  THIRD  PARAMETER    IS  OPTIMIZED  WITHIN   THE  MAIN 

C  PROGRAM  USING  AN   INCREMENTAL  OPTIMIZING  SEARCH  (i.e. 

C  FIRST  SEARCH   IN  TEN  PERCENT   INTERVALS  AND  THEN  TAKE 

C  MOST  PROMISING  TWENTY  PERCENT  REGION  AND  SEARCH  IN 

C  TWO  PERCENT   INTERVALS,  ETC.) 

C  OPTIMUM  OCCURS  WHEN  LOGARITHM  OF  LIKELIHOOD  FUNCTION 

C  IS  A  MAXIMUM 

C 

C  INPUT  DATA: 
C 

C  N--   NUMBER  OF   SAMPLES  -  ' 

C  NF--  NUMBER  OF  FAILURES 

C  I PROG-  FLAG  WHICH    INDICATES    IF  PROGRESSIVE 

C  CENSORING   IS   IN  EFFECT    (  1-- YES;  OTHER-NO) 

C  ICR-   NUMBER  OF  PARAMETERS  USED   IN  FIT 

C  (O-TWO5  <  >0-THREE) , 

C  T(I)-FAILURE  STRENGTHS  OF   SAMPLES  AS  CALCULATED  FOR  THE 

C  CENTRAL  PORTION  OF   THE  RING 

C  G( I )--CENSORED  VALUES  OF   STRENGTH    (READ    IN  FROM  SUBROUTINE 

C  IDC--FLAG  USED  TO  DETERMINE  STARTING  POINT  FOR 

C  OPTIMIZING  THIRD  PARAMETER    (AS  %  OF  1ST 

C  STRENGTH  VALUE) 

C  4-85>:;  OTHER-99/:: 

C  LOGICAL  UNIT  ASSIGNMENTS: 
C 

C  5-DATAFILE ( INPUT) 

C  9 CONSOLE 

C  6- DAT API LE (OUTPUT) 

C  /'-CONSOLE  (USED  BY  SUBROUTINE) 

C 


c***+>*-*  ***************** 

COMMON   / PAR AM /   ALPHA , BE T A , X  BAR  1 , S , N6 , X  NF 
D I  MENS I  ON  T ( 500 )  , TLOG ( 500 ) 

REAL  6 ( 1 50 ) , T 1 (1 50 ) , G 1 (1 50 ) , LSTAR (10) , G2 ( 1 50 ) , JR 

DATA  T/500*0. 0/ , TL0G/500*0. 0/ 

GFLAG-0 

READ (5,5)    N,NF  •  •  - 

5  FORMAT (213) 

READ (5,5)    I PROG, ICR 
AN  N 

DO  500   I    =    1 ,NF 
READ    (5,1)    T ( I ) 
1    FORMAT (F 15. 0)  ■  ' 

500  CONTINUE 

CALL  HAMLET ( N , NF , T , I PROG , G , GFLAG ) 
IF(ICR.EQ.O)    GOTO  223 

C 

C  BEGIN  OPTIMIZING  SEARCH  FOR  THIRD  PARAMETER 


WRITE (9, 15) 
1 5  FORMAT  (    '  I  NF'UT   I  DC    ( 3-  857. ,  4--95-/. )  ?  '  ) 

READ (9,*) IDC 

GFLAG=2 

DO  64   IV=1 ,NF 
64   Tl ( IV) =T (IV) 

NPG=M-NF 

DO  69  JV=1 ,NPG 
69  Gl ( JV) =G ( JV) 

XLHI--=99,  0 

IF    <IDC.EQ.3)  XLHI-95. 

IF    (IDC.EQ.4)  XLHI=85. 

XLLO^O. 0 

XL.STEP=--10. 
84  DO  93  JL=1 , 10 
93  LSTAR ( JL) =0. 0 

15=0 

C 

C  HERE   IS  THE  OPTIMIZATION  LOOP 

C 

DO  97    I R= 1 1 , 1 , -  1 

JR=XLHH-I5*XLSTEP 

15=15+1 

C= (Tl ( 1 ) *JR/100. ) 
IF(C.LT.O.O)  C=0.0 

C 

C  REEVALUATE  FAILURE  TIMES  AND  CENSORED  TIMES 

C 

DO  87   JQ=1 ,NF 
87   T ( JQ) =T1 ( JQ) -C 

DO  77   JP1=1 ,NP6 
77  G2  ( JPl )  =G1  (JPl )  -C 

C 

C  CALL   HAMLET   AGAIN  TO   COMPUTE   NEW    1ST   AND  2ND 

C  PARAMETERS 

C 

CALL   HAMLET ( N , NF , T , I PROG , G2 , GFLAG ) 

A=ALPHA 

B=BETA 

XNF=NF 

SUM 1=0. 0 

SUM2=0. 0 

DO  74  JZ=1 , XNF 

SUM 1 =SUM 1 +L06 ( T ( J  Z )  ) 
74  SUM2=SUM2+ (T (JZ) /BETA) ** ALPHA 

DO   76  J4=1,NP6 
76  SUM2=SUM2+ (G2 ( J4) /BETA ) **ALPHA 

C 

C  COMPUTE  LOGARITHM  OF  LIKELIHOOD  FUNCTION 

C 

154   LSTAR  (  15)  =XNF*LOG  (A)  -XNF*A*L06  (B)  +  (A-1  .  O )  *SUM  1 --SUM2 
WR ITE ( 6 , * ) LSTAR (15) , A , B , C , SUM 1 , SUM2 , NPG 
IF    ( 15. LE. 1 )    GOTO  97 
IF ( IDC. EQ. 1 )    GOTO  97 

I F    ( LSTAR (15). GT . LSTAR ( I 5- 1 ) )    GOTO  97 

C 

C  OPTIMUM   HAS   BEEN   LOCATED;    REDEFINE   SEARCH  INTERVAL 

C 

XNLO=JR 

XNHI  =  XNL0-i-2.  *ABS  (XLSTEP) 

IF    ( XNH I . GT. XLH I )  XNHI=XLHI 

IF    (XNLO.EQ. XNHI)    XNLO=XNH I+XLSTEP 

XLLO=XNLO 

XLHI=XNHI 

XLSTEP= ( XLLO-XLHI ) / 10. 


CHECK  FOR  PRECISION  OF    .0001  PERCENT 


97 


C1=ABS (LSTAR ( 15) -LSTAR (15- 
C2=ABS ( . 00000 1*LSTAR (15) ) 
IF    (C1.LE.C2)    SOTO  223 
GOTO  84 
CONTINUE 

OUTPUT  RESULTS 


1  )  ) 


223  WRITE (6, 6057) N,NF, BETA, ALPHA, C, XBARl ,S 

6057  FORMAT(///// '    NUMBER  OF  SPECIMENS  LOADED',  129// 
1'    NUMBER  OF  SPECIMENS  FAILED ', 129// 

WEIBULL   SCALE   PARAMETER ', F30. 5// 
WEIBULL   SHAPE   PARAMETER ', F30 . 5/ / 
6'    LOCATION  PARAMETER ' ,F30. 5// 
4'    LOBNORMAL  SCALE  PARAMETER ', F27 . 5/ / 
5'    LOBNORMAL  SHAPE   PARAMETER ', F27 . 5 ) 
WR I TE ( 6 , 6058 ) N , NF , ( T ( I ) , I = 1 , NF ) 

6058  FORMAT (213,1 000 ( / F 1 0 . O ) ) 
STOP 

END 

SUBROUT I NE  HAMLET ( N , NF , F , I  PROG , G , GFLAG ) 

COMMON   /PARAM/   ALPHA , BETA , XBAR 1 , S , NB , XNF 

REAL  F ( 150) ,6( 150) , X ( 150) 

REAL  MUU 

VUA  =1. 

X31  =9. 

ZC=0 

GIG  =1. 

IF ( IPROG. EQ. 1 )  GIG=0 

NFF  =  NF 

XNF  =  NF 

XN  =  N 

NG  =  N  -  NF 

NGG=NG 

I F ( NG . EQ . 0 )    GO  TO  511 0 
IF (GFLAG. NE. 0)    GOTO  5110 
DO  5 1 30   I   =   1 , NG 
IF (IPROG. EQ. 1 )    GOTO  5120 
G ( I )    =  F (NF) 
GOTO  5130 
READ (5,5121 ) 
FORMAT (FIO. 0) 
CONTINUE 
CONTINUE 


G  (  I  ) 


5120 
5121 
5130 
5 1 1 0 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  ORDER  FAILURE  AND  CENSORED  TIME  DATA  C 

C  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccc 

BAG  =  -10. 
GNG  =  NG 

IF    (GNG-.Ol)    1100,    1100,  1109 
1109      IF    (BIG  -   .01)    1106,    1106,  1100 
1106  NG1=NG-1 

DO    1112    1=1,  NGl 
11=1+1 

DO     1115  J=I 1 ,NG 

IF    (G(I)    -  G(J))    1115,    1115,  1118 
1118     GMID  =  G(J) 

B(J)    =  G(I) 

G(I)    =  GMID 
1115  CONTINUE 

WRITE(7.4887)  G(I) 


26. 


1 02 . 

134. 
135. 
136. 
141 . 
142. 


149. 


151 . 
152. 
153. 
154. 
1 55 . 
156. 
157. 
158. 
159. 
160. 
161 . 

162 

16 
164. 
165. 
166, 
167. 
168. 
1  A9. 


4887  FORMAK'O','      G(I)',F10,3)  170. 

1112  CONTINUE  171. 
WRITE (7,4887)  e(NG) 

1100     NFMl   =NF   -   1  172. 

DO  3    1=1 ,NFM1  173. 

11=1+1  174. 

DO   1   J=I1,NF  175. 

IF    (F<I)-F(J))    1,1,2  176. 

2  FMID=F(J)  177. 
F(J)=F(I)  178. 
F(I)=FMID  179. 

1            CONTINUE  180. 

X < I ) =AL0G10 (F ( I ) )  LOGGING 

3  CONTINUE  182. 
X(NF)    =  ALDG10(F (NF) )  183. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  1 84 . 

C  C  185. 

C          FIRST  TWO  ORDERED  FAILURE  TIMES  C  186. 

C  C  187. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  188. 

FNl    =  F  (1 )  189. 

FN2  =  F<2)  190. 

FNL  =  F<NF)  191. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  200, 

C  C  201. 

C           UNBIASED  ML  POINT  ESTIMATES  FOR  LOG  NORMAL  C  202. 

C            POP.    PARA.    ARE  CALCULATED   BELOW.  C  203. 

C  C  204. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  205 . 

XSUM=0.  SIGMA 

X2SUM=0.  XBAR 

DO  4   1=1 ,    NF  208. 

XSUM=XSUM+X ( I )  209. 

X2SUM=X2SUM+X ( I ) *X ( I )  210. 

4  CONTINUE  211. 
XN=N  212. 
XNF=NF  213. 
FNF  =  XNF  214. 
FN=  N  215. 
FN6  =  FN  -  FNF  216. 
XBAR=XSUM/XNF  217. 
XBARL  =  XBAR  218. 
S2=(X2SUM-XSUM*XBAR) / (XNF-1. )  219. 
S=SQRT(S2)  220. 
S=S+. 00001  221. 

cccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  C  255. 

C            UNBIASED  ESTIMATES  OF  LOG  NORMAL  MODEL   ARE          C  256. 

C           GIVEN  BELOW  FOR  CENSORED  CASE  C  257. 

C  C  258. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  259 . 

C            THIS  SECTION  FOR  HAMLET              CALLED    INSERT  C  260. 

IF  (FN  -  FNF  -  .01)  9100,  9100,  9103  261. 
9107  SIG  =.557 /ALPHA 

MUU  =  BETA  *   .56  **    (1.    /   ALPHA)  263. 

C            AAA  264. 

C            BBB  265. 

MUU  =  ALOGIO(MUU)  266. 

C  THE  BELOW  WORK   IS  COMMON  TO  HAMLET  AND  L0G7  268. 

8048      IF    (GIG  -  0.1)    8020,    8020,    8024  269. 

C  270. 

C              LOG  NORMAL  MLE   —   SINGLE   CENSORING  271. 

C  272. 

8024     EPS  = (ALOGIO (G ( 1 ) )    -  MUU)    /   SIG  273. 

ZE   =  CNORML      (EPS,    --150.,    0.0,    1.0)  274. 

FEZE  =  EXP (-(EPS   *EPS/2.))    *  0.3989423  275. 


IE  =  FEZE  / 


C 
C 

C 

8020 


3032 
8028 


4016 

401.2 

4008 
4024 

4020 


4C>04 
4005 
4000 
8040 


8044 

aose 

8050 


ZEE 


FMC3  * 


(  1  . 
ZE 


ZE) 


EASY  =  ZEE 


A I 
Bl 
CI 


ZEE  * 
ZEE  + 
EPS  * 


*  EPS 

<ZE  -EPS) 
EPS   *   A I 

(ZEE  +  BI) 


BO  TO  8028 


L06     NORMAL  MLE 


PROGRESSIVE  CENSORING 


ZEE  =  O. 
EASY  -  O. 
AI   =  0, 
BI   =  O. 
CI    =  0, 

DO  8032   I CE  =   1 ,  NG 
EPS  =    ( ALOG 1 0 ( 6 ( I CE ) )  - 
FEZE  =  EXP (-(EPS  *EPS/2. 


ZE  -  CNORML 
ZE  =  FEZE/ 
ZEE  =  ZEE  + 
EASY  =  EASY 
AP=ZE  *  (ZE 


ZE) 


BP 
CP 
AI 
BI 
CI 


ZE 
EPS 
A I  + 
BI    +  BP 
CI    -I-  CP 


(EP 
( 1 . 
ZE 

+  EPS  * 
-  EPS) 
EPS  *  AP 
■X-    (ZE  +  BP) 
AP 


IbO, 


MUU )  / 
)  )    *  O. 
,    0 .  0 , 


SIG 

3989423 
1  .  0) 


ZE 


CONTINUE 
TST 
TST2 
SR2 
FL  = 


/  SIG 


(  1 , 


GL 

PL 

QL 

RL 

DL 

EH= 

EK= 


(XBAR  -  MUU) 
=   TST   *  TST 
(S2   /    (SIG   *  SIG) > 
TST  +   ZEE   /  FNF 
TST2  +  SR2  -   1.    +  EASY/FNF 
AI/FNF)    /  SIG 
TST   +   BI /FNF) /SIG 


1 . /FNF) 


(  1 


-  \ 

-  (3 
PL  * 

(GL  * 
(FL  * 


RL 
QL 
QL 


(TST2  +  SR2) 
-QL   *  QL 

-  FL   *  RL)  / 

-  GL   *  PL)  / 


1 


+  CI /FNF)    /  SIG 


DL 
DL 


CONTINUE 

IF    (ABS(EK)    -  SIG/: 
EK   =  EK/2. 
GO   TO  4016 
IF    (ABS(EH)    -  0. 
EH   ==   EH/ 2. 
BO   TO  4008 
CONTINUE 
BAB  ~  -6. 
TEl   =  SIG  +EK 
TE2   =  MUU  +EH 
ZC  =   ZC  +  1. 
IF    (ZC-53,)  4000, 
CONT I NUE 
WRITE (7,4005) 
FORMAT ( 'SEE  LINE  : 
GO   TO  308 
CONTINUE 

IF    (TEl)    8040,  8040 
SIG  =  SIG/2. 
BAB  ==  6. 
GO  TO  8058 
SIG  =  TEl 

IF(TE2)    8050,  8050, 

MUU=MUU/2. 

BAB=6. 


02)    4008,    4008,    40 1: 


)    4020 ,    4020 ,  4024 


40(1)0 ,  4004 


29    IN  HAMLET ' ) 


8044 


3054 


276 
277 
278 
279 
280 
281 
282 
283 
284 


287 
288 
289 
290 
291 
292 
293 
294 
295 
296 
297 
298 
299 
300 
301 
302 
303 
304 
305 
306 
307 
308 
309 
310 
311 
312 
313 
314 
315 

317 
318 
319 
320 
321 

324 
325 
326 
327 
328 
329 


334 
335 
336 

338 
339 
340 
341 


8054 
8060 
8049 


8064 


C  THE 


C 
C 

c 
c 

9404 


C 
C 
C 
C 
C 

9 1 00 


GO  TO  8060 
liUU  ==  TE2 

IF    (BAB)    8049,    8048,  8048 
ERROR  =  ABS    (EH)+  ABS  (EK) 
AERR  =  0.0004 

IF  (ERROR  -  AERR)  8064,  8064,  8048 
XBAR  =  10.  **  MUU 
S=SIG 
ABOVE 
XBARL 
BNF  = 
ANF  = 


WORK    IS  COMMON  TO  HAMLET  AND  L0G7 
-  MUU 
XBAR 

2.    *  SIG 


IF    (BAG)    9111,    9111,  9404 


MLE  WE  I  BULL 
PARAMETERS 


PARAMETERS  TOO  LARGE — USE  DEFAULT 


ALPHA  =  0.557   /  S 
BEEL   =   XBARL   +  0.2506 
BETA  --=    10.  BEEL 
GO  TO  9111 


/  ALPHA 


133 


COMPUTE  NEW  FANG  PARAMETERS 


AL  =  0.557/S 

BEEL=XBAR  +  0.2506816/AL 

BEE=   10.    **  BEEL 

XBAR=    10.    **  XBARL 

WUM=0. 

WUMH=0. 

DO   133   1=1,  NF 
WUM=WUM  +  F(I) 
WUMH=  WUMH  +    1 . /F ( I ) 
CONTINUE 


R) 

(ESS/BNF  ~   1 . ) ) 
SECTION 


CALCS  SHAPE- 


ESS  =  WUM/XNF 
REC  =  WUMH/XNF 
R  =  l./REC 
BNF  =  SORT (ESS  * 
ANF  =  SORT    (2.  » 
CMLE  OF  WBL  PARAS  FOLLOW-FIRST 
ALO=  AL  »  VUA 
60  TO  12 

9103     ALO  =  F(l)    /    (F(2)    -  F(l)  +5) 
AL  =  0.557/S 

BEEL  =  XBAR  +  0.2506  /  AL 
BEE  =   10.    **  BEEL 
ALO  =  AL 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  C 
C  UNBIASED  ESTIMATES  OF  WEIBULL  MODEL  ARE  C 

C  GIVEN  BELOW  FOR  THE  COMPLETE  AND  CENSORED 

C  CASE  C 

C  C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
12  IF    (ALO*X (NF) -150. )  207,207,208 

207       SFA=  0. 

SFAX=  O. 
SFAX2  =  0. 
DO  5   1=1,  NF 

P  =  F(I)  ** 
SFA=  SFA  +  P 
SFAX=  SFAX  +  P 
SFAX2=  SFAX2  + 
5  CONTINUE 


LAST 


ALO 

* 
P 


X ( I ) *2 
*   X  ( I ) 


302585 
*   X  ( I  ) 


302585*2. 302585 


342. 

343. 

344. 

345. 

346. 

347. 

348. 

349. 

350. 

351 . 

352. 

353. 

354. 

355. 

356. 

357. 

358. 

359. 

360. 

361 . 

362. 

363. 

364. 

365. 

366. 

367. 

368. 

369. 
BE 

371  . 
NWFGLD 

374 . 

375. 
376, 
377. 
378. 
379. 
380. 
BNF 
382. 

SECTION  FOR  NEW  FANG 
384. 
388. 
389. 
390 . 
391 . 
392. 
393. 
394. 
395. 
396. 
397. 
398. 
399. 
400 . 
401 . 
402. 
403. 
404. 
405. 
406. 
407. 
408. 
409. 
410. 


SGA=0.  411. 

SBAX=0.  412. 

SGAX2=0.  413. 

IF    (GNG   -    .01)    1130,    1130,    1133  414. 

1133      IF    (GNG  -  7.01)    1136,    1136,    1139  415. 

1139      IF    (GIG   -    .01)    1136,    1136,    1142  ,  416. 

1142     P==  G(l)    **  ALO  417. 

GLOGG  =   ALDG    (G(l))  418. 

SGA=   GNG   *   P  419. 

SGAX=  SGA  *  GLOGG  420. 

SGAX2=  SGAX   >   GLOGG  421. 

GO  TO   1130  .  422. 

1136     DO   1148   1=1,    NG  423. 

P=  G(I)    **ALO  424. 

GLOGG=  ALOG    (G(I))  .  425. 

SGA=   SGA  +   P  426. 

SGAX-  SGAX   +  GLOGG  *  P  427. 

SGAX2=  SGAX2  +  P  *  GLOGG  *  GLOGG  428. 

1148     CONTINUE  429. 

1130     SFA=   SFA   +   SGA  430. 

SFAX   =  SFAX   +  SGAX  v  431. 

SFAX2   =   SFAX2  +   SGAX2  432. 

FALO=  SFAX/SFA  -    1 . /ALO  -   XBARL*2 . 302585  433. 

SFA2-2*AL0G (SFA)  434 

DALO 1 = ALOG ( SFA  X  2 ) -ALOG ( SFA ) 

DAL02=2*AL0B (SFAX ) -SFA2 

DFAL0= (10. **DAL01 ) - (10. **DAL02) +1 . / (AL0*AL0) 

"32  IF    (ABS(DFALO)    -   .00001)    30,    30,    31  ^  436. 

30  DFALO=  10.  *  DFALO  437. 
GO  TO  32  438. 

31  CORR=  -FALO/ DFALO  439. 
IF    (ABS(CORR)         .0001)    8,    8,    9  440, 

9  TRY=  ALO  +  CORR  441. 
IF    (TRY)    10,    10,    11  442. 

10  ALO=  ALO/ 2.  443. 
GO   TO    12  444. 

11  ALO=  TRY  445. 
GO  TO   12  446. 

S  ALPHA=  ALO  +  CORR  447. 

BETA  =    (SFA/XNF)    **    (1. /ALPHA)  448. 

GO   TO  223  449. 

208  CONTINUE  450. 
IF    (FNG   -   0.1)    209,    210,    210  451. 

C  452. 

C  STATEMENT   12  POSITIVE — WE I  BULL  PARAMETERS  ABORTED  453. 

C  454. 

210  WRITE (7, 211) 

211  FORMAT ('SEE  LINE  453    IN  HAMLET') 
GO   TO  308 

209  CONTINUE  .  457. 
ALPHA==AL  458. 
BETA---BEE  ^  459. 
BAC   =  8.                                                                V  460. 

223       CONTINUE  461. 

I F ( FNG  -    .01)    9111,9111,91 07 

C  463. 

C  DETERMINE   IF  COMPLETE  SAMPLE:      IF  NOT  FIND  464. 

C  CENSORED  LOG  NORMAL   PARAMETERS  465. 

C  466. 

9111      CONTINUE  467. 

ALPHA  =  ALPHA/ VUA  468. 

NL=NFF  469. 

NW=NFF  470. 

XBARl    =   ALOG(XBAR)  488.1 

S   =   S*AL0G(10.)  488.2 

308       CONTINUE  504- 


RETURN 

END  505 

FUNCTION  CNORML.    (  XH  ,  XL. ,  XM  ,  XS )  506 

C            S00887        CNORML        F"'S-497          CHIANG   E.    C.    H.        661004   6600  507 

P(X)=.5*(1.-(1. / 141 12821*X+. 0B864027*X**2  508 

1-I-.  02743349*X**3-.  00039446* X**4-i- .  00328975*X**5  509 

2)**8))  510 

C                              IF   NM=1,X1    AND   X2  SHOULD  BE  STANDARDIZED  TO  N(0,1)  511 

Xl= <XH-XN) /XS  512 

X2=(XL-Xri) /XS  513 

Z1=X1/1 . 414213567  514 

XI   =   Zl  515 

Z2=X2/1 . 414213567  516 

IF (Z1*Z2) 1 ,2,3  517 

C            WHEN   Zl    AND   Z2  HAVE  DIFFERENT   SIGN  518 

1  Z2=ABS(Z2)  519 
CN0F^ML=P(Z2)+P(Z1)  520 
GO  TO    100  521 

C            TO  FIND  WHETHER   Zl    OR   Z2    IS  0  522 

2  CN0RML=P(ABS(Z1+Z2) )  523 
GO   TO    100  524 

C            WHEN   Zl    AND   Z2  HAVE  THE  SAME  SIGN  525 

3  Z2=ABS(Z2)  526 
Z1=-ABS(Z1)  527 
CN0RML-ABS(P(Z2)-P(Z1) )  528 

100        RETURN  529 

END  530 

t-BEND 
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